# DPIV_ALGORITHM
from scipy.interpolate import griddata
from scipy.io import savemat
import matplotlib.pyplot as plt
import numpy as np
import cv2
import os
import re
import shutil
import zipfile
import urllib.request
from dpivsoft import DPIV
# Image formats accepted by the stereo helpers.
IMAGE_EXTS = ('.tif', '.tiff', '.png', '.bmp', '.jpg', '.jpeg')
# 4th International PIV Challenge (2014), Case E: time-resolved stereo-PIV
# vortex ring (two cameras, 1 and 3). ~197 MB, hosted on pivchallenge.org.
PIVCHALLENGE_E_URL = (
"https://www.pivchallenge.org/pub14/4th_PIV-Challenge_Case_E.zip")
[docs]
def rgb2gray(item):
"""Return *item* as a 2-D float grayscale array.
Accepts either a file path (any format ``cv2`` can read) or an
already-loaded ndarray, so callers can mix disk files and in-memory
back-projected frames interchangeably.
"""
if isinstance(item, np.ndarray):
img = item
else:
img = cv2.imread(item, cv2.IMREAD_UNCHANGED)
if img is None:
raise FileNotFoundError(item)
if img.ndim == 3:
img = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)
return img.astype(float)
[docs]
def stereo_calibration(final_matrix_l, final_matrix_r, org_width, org_height,
mode='reduced'):
"""Calibrate a stereo camera pair using the Soloff mapping.
``mode`` selects the polynomial:
* ``'reduced'`` (default) — the 9-term, linear-in-z model
(basis [x², y², xy, x, y, 1, z, x·z, y·z]); needs >=2 distinct z planes.
* ``'full'`` — the classic 19-term Soloff (Soloff et al. 1997): cubic
in-plane at z⁰ (the cubic terms model lens distortion), quadratic at z¹,
linear at z². Needs >=3 distinct z planes (the highest z-power) and
enough in-plane points per plane. Falls back to NaN maps (Soloff
unavailable) if the plane condition is not met.
For each camera, fits the chosen polynomial in both directions by exact
least squares:
forward (pixel -> world): X,Y = poly(x_px, y_px, Z)
map (world -> pixel): x_px,y_px = poly(X, Y, Z)
The forward map's Z=0 part builds the world-coordinate grid; the world->pixel
`map` is the Soloff mapping used by soloff_velocity (full local Jacobian) and
for back-projection (its Z=0 part). A dual-plane calibration (>=2 distinct Z)
is required to constrain the parallax (Z and cross-terms).
Returns
-------
fwd_a_l, fwd_b_l, fwd_a_r, fwd_b_r : ndarray, shape (9,) or (19,)
Forward (pixel→world) coefficients for X and Y, per camera (9 in reduced
mode, 19 in full mode). Used (Z=0 part) to compute the world grid.
map_a_l, map_b_l, map_a_r, map_b_r : ndarray, shape (9,) or (19,)
Soloff world->pixel mapping coefficients for x and y, per camera; used by
soloff_velocity (local Jacobian) and back-projection (Z=0 part). NaN-filled
if the calibration lacks enough distinct Z planes (parallax unconstrained).
results : str
Human-readable summary of fit RMS and Soloff-map readiness.
"""
full = (mode == 'full')
min_planes = 3 if full else 2
_fit = _fit_poly19 if full else _fit_poly9
_design = _poly19_design if full else _poly9_design
n_terms = 19 if full else 9
def _fit_camera(x_px, y_px, X_w, Y_w, Z_w):
"""Fit both Soloff maps for one camera by least squares.
Returns (fwd_a, fwd_b, map_a, map_b), each an ``n_terms``-vector:
* fwd : pixel -> world (X,Y = fwd(x_px, y_px, Z)); its Z=0 part
fwd[:6] is the in-plane map used to build the world grid.
* map : world -> pixel (the Soloff map; x_px,y_px = map(X,Y,Z)), used
by soloff_velocity (full local Jacobian) and for back-projection
(its Z=0 part map[:6]).
Needs >=``min_planes`` distinct Z planes to constrain the z (and, in full
mode, z²) terms; otherwise returns NaN maps.
"""
if len(np.unique(np.round(Z_w, 9))) < min_planes:
nan = lambda: np.full(n_terms, np.nan)
return nan(), nan(), nan(), nan()
fwd_a = _fit(x_px, y_px, Z_w, X_w) # pixel (+Z) -> world X
fwd_b = _fit(x_px, y_px, Z_w, Y_w) # pixel (+Z) -> world Y
map_a = _fit(X_w, Y_w, Z_w, x_px) # world (+Z) -> pixel x
map_b = _fit(X_w, Y_w, Z_w, y_px) # world (+Z) -> pixel y
return fwd_a, fwd_b, map_a, map_b
x_l, y_l = final_matrix_l[:, 0], final_matrix_l[:, 1]
X_l, Y_l, Z_l = final_matrix_l[:, 2], final_matrix_l[:, 3], final_matrix_l[:, 4]
x_r, y_r = final_matrix_r[:, 0], final_matrix_r[:, 1]
X_r, Y_r, Z_r = final_matrix_r[:, 2], final_matrix_r[:, 3], final_matrix_r[:, 4]
print("Fitting left camera polynomial...")
fwd_a_l, fwd_b_l, map_a_l, map_b_l = _fit_camera(x_l, y_l, X_l, Y_l, Z_l)
print("Fitting right camera polynomial...")
fwd_a_r, fwd_b_r, map_a_r, map_b_r = _fit_camera(x_r, y_r, X_r, Y_r, Z_r)
def _rms_fwd(a, b, x_px, y_px, Z_w, X_w, Y_w):
"""Forward-fit residual (pixel->world) for one camera, RMS in mm."""
C = _design(x_px, y_px, Z_w)
return float(np.sqrt(np.mean((C.dot(a) - X_w)**2 + (C.dot(b) - Y_w)**2)))
rms_l = _rms_fwd(fwd_a_l, fwd_b_l, x_l, y_l, Z_l, X_l, Y_l)
rms_r = _rms_fwd(fwd_a_r, fwd_b_r, x_r, y_r, Z_r, X_r, Y_r)
# Soloff map reprojection residual (world->pixel, in pixels) and validity:
# this is what soloff_velocity uses. NaN maps mean the calibration had a
# single Z plane (no parallax) and Soloff is not possible.
def _rms_map(map_a, map_b, X_w, Y_w, Z_w, x_px, y_px):
"""Soloff map reprojection residual (world->pixel), RMS in px;
NaN when the map is NaN (single Z plane, no parallax)."""
if np.any(np.isnan(map_a)) or np.any(np.isnan(map_b)):
return float('nan')
C = _design(X_w, Y_w, Z_w)
return float(np.sqrt(np.mean(
(C.dot(map_a) - x_px)**2 + (C.dot(map_b) - y_px)**2)))
map_rms_l = _rms_map(map_a_l, map_b_l, X_l, Y_l, Z_l, x_l, y_l)
map_rms_r = _rms_map(map_a_r, map_b_r, X_r, Y_r, Z_r, x_r, y_r)
soloff_ok = not (np.isnan(map_rms_l) or np.isnan(map_rms_r))
# Diagnostic plot: FOV box corners in world coords, at the Z=0 plane.
m = 48
cx = np.array([m, org_width-m, org_width-m, m, m], dtype=float)
cy = np.array([m, m, org_height-m, org_height-m, m], dtype=float)
Xc_l, Yc_l = _poly_eval(fwd_a_l, cx, cy, 0.0), _poly_eval(fwd_b_l, cx, cy, 0.0)
Xc_r, Yc_r = _poly_eval(fwd_a_r, cx, cy, 0.0), _poly_eval(fwd_b_r, cx, cy, 0.0)
# Reprojection check: project the measured calibration pixels (at their own
# Z) to world coords through each camera's forward map. How tightly these
# land on the true calibration points -- and on each other -- is the visual
# confirmation of the fit; their spread is the forward RMS in the title.
X_proj_l = _poly_eval(fwd_a_l, x_l, y_l, Z_l)
Y_proj_l = _poly_eval(fwd_b_l, x_l, y_l, Z_l)
X_proj_r = _poly_eval(fwd_a_r, x_r, y_r, Z_r)
Y_proj_r = _poly_eval(fwd_b_r, x_r, y_r, Z_r)
fig, ax = plt.subplots()
ax.plot(Xc_l, Yc_l, '*-', label='Left camera field of view')
ax.plot(Xc_r, Yc_r, 'o-', label='Right camera field of view')
ax.plot(np.concatenate([X_l, X_r]),
np.concatenate([Y_l, Y_r]), '.k', label='Calibration points (true)')
ax.plot(X_proj_l, Y_proj_l, 'o', mfc='none', color='tab:blue',
label='Left camera projected points')
ax.plot(X_proj_r, Y_proj_r, 'x', color='tab:orange',
label='Right camera projected points')
ax.plot(0, 0, '+k', markersize=25, label='Origin')
ax.set_xlabel('X (mm)')
ax.set_ylabel('Y (mm)')
ax.set_aspect('equal')
ax.set_title(
f'Polynomial calibration — RMS: left {rms_l:.3f} mm right {rms_r:.3f} mm')
ax.legend()
plt.show()
model_name = "Full (19-term, cubic)" if full else "Reduced (9-term)"
if soloff_ok:
soloff_line = (
f" Soloff maps: OK — reprojection RMS "
f"L {map_rms_l:.3f} px / R {map_rms_r:.3f} px (Soloff-ready)"
)
else:
soloff_line = (
f" Soloff maps: NOT available — needs >={min_planes} distinct Z "
f"planes for the {model_name} model.\n"
" Use a dual-plane TSI/LaVision target (or merge multi-depth CPTs)."
)
results = (
f" Model: {model_name} Soloff\n"
f" Forward fit RMS — Left: {rms_l:.3f} mm Right: {rms_r:.3f} mm\n"
+ soloff_line
)
return (fwd_a_l, fwd_b_l, fwd_a_r, fwd_b_r,
map_a_l, map_b_l, map_a_r, map_b_r, results)
[docs]
def _poly9_design(px, py, pz):
"""Reduced (9-term) Soloff design matrix.
Basis [p², q², pq, p, q, 1, z, pz, qz] — quadratic in-plane at z⁰, linear
in-plane at z¹. KEEP this order (on-disk files and the hand-derived
fold-back in _fit_poly9 depend on it).
"""
n = len(px)
C = np.zeros((n, 9))
C[:, 0] = px ** 2
C[:, 1] = py ** 2
C[:, 2] = px * py
C[:, 3] = px
C[:, 4] = py
C[:, 5] = 1.0
C[:, 6] = pz
C[:, 7] = px * pz
C[:, 8] = py * pz
return C
[docs]
def _poly19_design(px, py, pz):
"""Full (19-term) Soloff design matrix (Soloff et al. 1997).
Basis, in coefficient order: cubic in-plane at z⁰ (the cubic terms model
lens distortion), quadratic in-plane at z¹, linear in-plane at z²:
[p³, p²q, pq², q³, p², pq, q², p, q, 1, (z⁰, 10)
p²z, pqz, q²z, pz, qz, z, (z¹, 6)
pz², qz², z²] (z², 3)
"""
n = len(px)
C = np.zeros((n, 19))
# z⁰: cubic in-plane
C[:, 0] = px ** 3
C[:, 1] = px ** 2 * py
C[:, 2] = px * py ** 2
C[:, 3] = py ** 3
C[:, 4] = px ** 2
C[:, 5] = px * py
C[:, 6] = py ** 2
C[:, 7] = px
C[:, 8] = py
C[:, 9] = 1.0
# z¹: quadratic in-plane
C[:, 10] = px ** 2 * pz
C[:, 11] = px * py * pz
C[:, 12] = py ** 2 * pz
C[:, 13] = px * pz
C[:, 14] = py * pz
C[:, 15] = pz
# z²: linear in-plane
C[:, 16] = px * pz ** 2
C[:, 17] = py * pz ** 2
C[:, 18] = pz ** 2
return C
[docs]
def _fit_poly9(px, py, pz, b):
"""Least-squares fit of the 9-term Soloff polynomial, well-conditioned.
Fits ``b ≈ poly9(px, py, pz)`` but solves in *normalised* coordinates (each
of px, py, pz shifted to zero mean and unit std), so the design-matrix
columns are all O(1) instead of spanning ~7 orders of magnitude (px² in
pixels ~1e7, the z column ~1e0). That drops cond(C) from ~1e7 to ~1, so the
weakly-excited parallax coefficients (z, p·z, q·z) — the ones stereo depends
on — are resolved at full precision.
The fitted coefficients are folded back into the *raw* 9-term basis, so the
returned 9-vector is in the exact same convention as a direct raw fit (the
basis is closed under per-axis affine normalisation): ``_poly9_eval`` /
``_poly9_jacobian``, every caller, and the on-disk format stay unchanged.
"""
px = np.asarray(px, dtype=float)
py = np.asarray(py, dtype=float)
pz = np.asarray(pz, dtype=float)
mx, my, mz = px.mean(), py.mean(), pz.mean()
sx = px.std() or 1.0
sy = py.std() or 1.0
sz = pz.std() or 1.0
C = _poly9_design((px - mx) / sx, (py - my) / sy, (pz - mz) / sz)
c = np.linalg.lstsq(C, b, rcond=None)[0] # coeffs in normalised basis
# Substitute p' = (p - m)/s (and q, z likewise) into the normalised
# polynomial and collect terms — the 9-term basis maps back onto itself.
sx2, sy2 = sx * sx, sy * sy
raw = np.empty(9)
raw[0] = c[0] / sx2
raw[1] = c[1] / sy2
raw[2] = c[2] / (sx * sy)
raw[3] = -2*c[0]*mx/sx2 - c[2]*my/(sx*sy) + c[3]/sx - c[7]*mz/(sx*sz)
raw[4] = -2*c[1]*my/sy2 - c[2]*mx/(sx*sy) + c[4]/sy - c[8]*mz/(sy*sz)
raw[5] = (c[0]*mx*mx/sx2 + c[1]*my*my/sy2 + c[2]*mx*my/(sx*sy)
- c[3]*mx/sx - c[4]*my/sy + c[5] - c[6]*mz/sz
+ c[7]*mx*mz/(sx*sz) + c[8]*my*mz/(sy*sz))
raw[6] = c[6]/sz - c[7]*mx/(sx*sz) - c[8]*my/(sy*sz)
raw[7] = c[7]/(sx*sz)
raw[8] = c[8]/(sy*sz)
return raw
[docs]
def _fit_poly19(px, py, pz, b):
"""Least-squares fit of the full 19-term Soloff polynomial.
Cubic in-plane at z⁰ (lens distortion), quadratic at z¹, linear at z²
(Soloff et al. 1997). Needs >=3 distinct z planes (the highest z-power) and
enough in-plane points per plane for the cubic terms.
Conditioning is handled by *column equilibration* rather than the
coordinate-normalisation fold-back used by ``_fit_poly9``: each design
column is scaled to unit RMS before the solve and the coefficients are
divided back afterwards (a similarity transform that leaves the raw
least-squares solution unchanged but drops cond(C) by many orders).
Equilibration is basis-agnostic, so it extends to all 19 terms without
re-deriving the fold-back algebra.
"""
px = np.asarray(px, dtype=float)
py = np.asarray(py, dtype=float)
pz = np.asarray(pz, dtype=float)
C = _poly19_design(px, py, pz)
scale = np.sqrt(np.mean(C**2, axis=0))
scale[scale == 0] = 1.0
c = np.linalg.lstsq(C / scale, b, rcond=None)[0]
return c / scale
[docs]
def _is_full_poly(c):
"""True for a full (19-term) Soloff vector, False for reduced (9-term).
Calibrations are only ever 9 (reduced) or 19 (full Soloff); any other length
is a malformed/corrupt map, so raise rather than silently mis-modelling.
Used by the ``_poly_eval`` / ``_poly_jacobian`` selectors to dispatch.
"""
n = len(c)
if n == 19:
return True
if n == 9:
return False
raise ValueError(
f"Soloff coefficient vector must have length 9 (reduced) or 19 (full), "
f"got {n}.")
# --- Reduced (9-term) model -------------------------------------------------
[docs]
def _poly9_eval(c, X, Y, Z):
"""Evaluate the reduced 9-term Soloff polynomial.
Basis [X², Y², XY, X, Y, 1, Z, X·Z, Y·Z]. X, Y, Z may be arrays.
"""
return (c[0]*X**2 + c[1]*Y**2 + c[2]*X*Y + c[3]*X + c[4]*Y + c[5]
+ c[6]*Z + c[7]*X*Z + c[8]*Y*Z)
[docs]
def _poly9_jacobian(c, X, Y, Z=0.0):
"""Partials (dF/dX, dF/dY, dF/dZ) of the reduced 9-term polynomial at Z.
Linear in Z, so the parallax slope dF/dZ = c6 + c7·X + c8·Y does NOT depend
on Z; only the in-plane derivatives pick up the small c7·Z / c8·Z terms.
"""
dFdX = 2*c[0]*X + c[2]*Y + c[3] + c[7]*Z
dFdY = 2*c[1]*Y + c[2]*X + c[4] + c[8]*Z
dFdZ = c[6] + c[7]*X + c[8]*Y
return dFdX, dFdY, dFdZ
# --- Full (19-term) Soloff model --------------------------------------------
[docs]
def _poly19_eval(c, X, Y, Z):
"""Evaluate the full 19-term Soloff polynomial.
Coefficient order matches ``_poly19_design``: cubic in-plane at z⁰,
quadratic at z¹, linear at z². X, Y, Z may be arrays.
"""
return (
# z⁰: cubic in-plane
c[0]*X**3 + c[1]*X**2*Y + c[2]*X*Y**2 + c[3]*Y**3
+ c[4]*X**2 + c[5]*X*Y + c[6]*Y**2 + c[7]*X + c[8]*Y + c[9]
# z¹: quadratic in-plane
+ c[10]*X**2*Z + c[11]*X*Y*Z + c[12]*Y**2*Z
+ c[13]*X*Z + c[14]*Y*Z + c[15]*Z
# z²: linear in-plane
+ c[16]*X*Z**2 + c[17]*Y*Z**2 + c[18]*Z**2)
[docs]
def _poly19_jacobian(c, X, Y, Z=0.0):
"""Partials (dF/dX, dF/dY, dF/dZ) of the full 19-term polynomial at Z.
Unlike the reduced model, the parallax slope dF/dZ is Z-dependent (the z²
curvature the extra planes buy) and the in-plane derivatives carry the
cubic-distortion terms.
"""
dFdX = (3*c[0]*X**2 + 2*c[1]*X*Y + c[2]*Y**2 + 2*c[4]*X + c[5]*Y + c[7]
+ 2*c[10]*X*Z + c[11]*Y*Z + c[13]*Z + c[16]*Z**2)
dFdY = (c[1]*X**2 + 2*c[2]*X*Y + 3*c[3]*Y**2 + c[5]*X + 2*c[6]*Y + c[8]
+ c[11]*X*Z + 2*c[12]*Y*Z + c[14]*Z + c[17]*Z**2)
dFdZ = (c[10]*X**2 + c[11]*X*Y + c[12]*Y**2 + c[13]*X + c[14]*Y + c[15]
+ 2*c[16]*X*Z + 2*c[17]*Y*Z + 2*c[18]*Z)
return dFdX, dFdY, dFdZ
# --- Length-dispatch selectors (used by every caller) -----------------------
[docs]
def _poly_eval(c, X, Y, Z):
"""Evaluate the Soloff map, picking the model by coefficient length."""
if _is_full_poly(c):
return _poly19_eval(c, X, Y, Z)
return _poly9_eval(c, X, Y, Z)
[docs]
def _poly_jacobian(c, X, Y, Z=0.0):
"""Soloff map Jacobian, picking the model by coefficient length."""
if _is_full_poly(c):
return _poly19_jacobian(c, X, Y, Z)
return _poly9_jacobian(c, X, Y, Z)
[docs]
def _stereo_recon_grid(calibration, img_shape, piv_shape):
"""Common world grid at PIV vector resolution over the cameras' overlap.
Shared by ``soloff_velocity`` and ``disparity_to_delta_z`` so both work
on identical grid points. ``make_world_grid`` returns the overlap at full
image resolution (for dewarping); reconstructing there would upsample a few
thousand PIV vectors onto millions of points -- needlessly huge and
physically meaningless -- so the grid is sampled at the input PIV shape.
Returns (Xg, Yg): 1-D world axes (mm) meshgrid.
"""
wx, wy = make_world_grid(calibration, img_shape)
if len(piv_shape) == 2:
ny_piv, nx_piv = piv_shape
else:
n = max(2, int(round(np.sqrt(int(np.prod(piv_shape))))))
ny_piv = nx_piv = n
world_x = np.linspace(wx[0], wx[-1], nx_piv)
world_y = np.linspace(wy[0], wy[-1], ny_piv)
Xg, Yg = np.meshgrid(world_x, world_y)
return Xg, Yg
[docs]
def soloff_velocity(x_l, y_l, u_l, v_l, x_r, y_r, u_r, v_r, calibration,
img_shape, delta_z=None):
"""Classic Soloff 3-component reconstruction from raw-image PIV fields.
At every point of a common world grid it builds each camera's *local*
world->pixel mapping Jacobian and solves
the over-determined linear system
[Δx_l, Δy_l, Δx_r, Δy_r]ᵀ = M · [ΔX, ΔY, ΔZ]ᵀ (M is 4×3)
by least squares for the world displacement (U, V, W). The Jacobian (hence
the parallax) varies across the field, so the reconstruction is locally
correct everywhere.
Parameters
----------
x_l, y_l : ndarray — left-camera PIV grid in *raw image* pixels.
u_l, v_l : ndarray — left-camera in-plane displacement (raw pixels).
x_r, y_r, u_r, v_r : ndarray — same for the right camera.
calibration : Parameters.Stereo_Calibration
Must carry the Soloff maps ``map_a_l/map_b_l/map_a_r/map_b_r`` (produced
by a dual-plane calibration).
img_shape : tuple (height, width)
Pixel dimensions of the raw camera images (for the world grid extent).
delta_z : None, float, or ndarray, optional
Out-of-plane offset Z (mm) at which to evaluate the Soloff map and its
Jacobian, instead of the calibration plane Z=0. This is the disparity
/ Δz self-calibration correction (Wieneke 2005): it accounts for the
laser sheet sitting at a (possibly spatially varying) offset from the
calibration plane. ``None`` (default) reproduces the uncorrected Z=0
reconstruction. A scalar applies a uniform offset; an array must match
the input PIV grid shape (``np.shape(x_l)``) and gives a per-point Δz —
build it from a measured disparity field with ``disparity_to_delta_z``.
Because the map already carries the Z dependence.
Returns
-------
x, y : ndarray — common world grid (mm).
u, v : ndarray — in-plane world displacement.
u_z : ndarray — out-of-plane (w) displacement. Front-Back convention:
negative = toward left camera, positive = toward right camera.
"""
maps = {'map_a_l': getattr(calibration, 'map_a_l', None),
'map_b_l': getattr(calibration, 'map_b_l', None),
'map_a_r': getattr(calibration, 'map_a_r', None),
'map_b_r': getattr(calibration, 'map_b_r', None)}
for name, c in maps.items():
if c is None or np.any(np.isnan(np.asarray(c, dtype=float))):
raise ValueError(
"Soloff reconstruction requires dual-plane calibration maps "
f"('{name}' missing or NaN).")
map_a_l, map_b_l = maps['map_a_l'], maps['map_b_l']
map_a_r, map_b_r = maps['map_a_r'], maps['map_b_r']
# Common world grid over the cameras' overlapping field of view, sampled
# at the PIV VECTOR resolution (see _stereo_recon_grid).
Xg, Yg = _stereo_recon_grid(
calibration, img_shape, np.shape(x_l))
shape = Xg.shape
Xf, Yf = Xg.ravel(), Yg.ravel()
# Out-of-plane evaluation plane: Z=0 (calibration plane) unless a disparity
# Δz correction is supplied. The map carries the Z dependence, so this
# shifts where each world point images / its parallax -- no image re-warp.
if delta_z is None:
Zf = 0.0
else:
dz = np.asarray(delta_z, dtype=float)
if dz.ndim == 0:
Zf = float(dz)
elif dz.size == Xf.size:
Zf = dz.ravel()
else:
raise ValueError(
f"delta_z shape {dz.shape} does not match the reconstruction "
f"grid {shape} (= input PIV grid shape np.shape(x_l)).")
def _cam_disp(map_a, map_b, xg, yg, u, v):
"""Raw-image displacement of each world grid point for one camera."""
x_img = _poly_eval(map_a, Xf, Yf, Zf)
y_img = _poly_eval(map_b, Xf, Yf, Zf)
src = (xg.ravel(), yg.ravel())
dpx = griddata(src, u.ravel(), (x_img, y_img), method='linear')
dpy = griddata(src, v.ravel(), (x_img, y_img), method='linear')
return dpx, dpy
dpx_l, dpy_l = _cam_disp(map_a_l, map_b_l, x_l, y_l, u_l, v_l)
dpx_r, dpy_r = _cam_disp(map_a_r, map_b_r, x_r, y_r, u_r, v_r)
# Per-camera mapping Jacobian rows (each ∂(x or y)/∂(X,Y,Z)) at every point,
# evaluated at the (possibly disparity-corrected) plane Z=Zf.
jxl = _poly_jacobian(map_a_l, Xf, Yf, Zf)
jyl = _poly_jacobian(map_b_l, Xf, Yf, Zf)
jxr = _poly_jacobian(map_a_r, Xf, Yf, Zf)
jyr = _poly_jacobian(map_b_r, Xf, Yf, Zf)
# M: (N, 4, 3) stacked [Δx_l; Δy_l; Δx_r; Δy_r] = M·[ΔX, ΔY, ΔZ]
M = np.stack([np.column_stack(jxl), np.column_stack(jyl),
np.column_stack(jxr), np.column_stack(jyr)], axis=1)
d = np.stack([dpx_l, dpy_l, dpx_r, dpy_r], axis=1) # (N, 4)
# Solve normal equations per point; mask points where any camera had o
# valid interpolated displacement (outside its field of view).
valid = np.all(np.isfinite(d), axis=1)
d_filled = np.where(np.isfinite(d), d, 0.0)
MtM = np.einsum('nij,nik->njk', M, M) # (N, 3, 3)
Mtd = np.einsum('nij,ni->nj', M, d_filled) # (N, 3)
sol = np.linalg.solve(MtM, Mtd) # (N, 3)
sol[~valid] = np.nan
u = sol[:, 0].reshape(shape)
v = sol[:, 1].reshape(shape)
u_z = sol[:, 2].reshape(shape)
# Front-Back w convention: negative = toward left camera, positive = toward
# right. With the corrected get_real_TSI Front-Back Z assignment the
# calibration frame already has +Z toward the right camera, so the raw ΔZ
# from the least-squares solve is already in this convention (no flip).
return Xg, Yg, u, v, u_z
[docs]
def save(x, y, u, v, u_z, filename, option='dpivsoft', Matlab=False, param=False):
"""
save flow field to a file. Option indicates the saving
format.
dpivsof: save in python .npz file using the original
formating of dpivsoft in matlab
openpiv: save the field in an ascii file compatible
with openpiv
"""
from dpivsoft.Classes import Parameters
u = u / Parameters.delta_t
v = v / Parameters.delta_t
u_z = u_z / Parameters.delta_t
if Matlab:
mdic = {"x": x*1.0, "y": y*1.0, "u": u*1.0, "v": v*1.0, "u_z": u_z*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),
"mask": float(Parameters.mask)}
savemat(filename + '.mat', mdic)
if option == 'dpivsoft':
if param:
np.savez(filename, x=x, y=y, u=u, v=v, u_z=u_z,
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, u_z=u_z,
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, u_z, np.ones_like(x, dtype=int)]])
np.savetxt(
filename,
out.T,
fmt=fmt,
delimiter=delimiter,
header="x"
+ delimiter
+ "y"
+ delimiter
+ "u"
+ delimiter
+ "v"
+ delimiter
+ "u_z"
+ delimiter
+ "mask",
)
else:
raise ValueError(f"Unknown saving option: '{option}'")
[docs]
def make_world_grid(calibration, img_shape, scale_px_per_mm=None, margin_mm=0.0):
"""Compute a world-coordinate output grid covering both cameras' fields of view.
Maps the image corners of each camera to world coordinates via the forward
polynomial mapping, then builds a regular grid over the intersection bounding box.
Parameters
----------
calibration : Parameters.Stereo_Calibration
img_shape : tuple (height, width)
Pixel dimensions of the raw camera images.
scale_px_per_mm : float or None
Output resolution (pixels per mm). If None, estimated from the
average image-width-to-world-X-extent ratio of both cameras.
margin_mm : float
Inward margin subtracted from each side of the overlap region (mm).
Returns
-------
world_x : ndarray, shape (nx,) — X coordinates of grid columns (mm).
world_y : ndarray, shape (ny,) — Y coordinates of grid rows (mm),
running from y_max (row 0, physically high) to y_min (last row).
"""
h, w = img_shape
corners_x = np.array([0.0, w, w, 0.0])
corners_y = np.array([0.0, 0.0, h, h])
# Forward map (pixel->world) at the Z=0 plane; the selector handles the
# 9- or 19-term in-plane part (the cubic terms matter for the 19-term map).
Xl = _poly_eval(calibration.fwd_a_l, corners_x, corners_y, 0.0)
Yl = _poly_eval(calibration.fwd_b_l, corners_x, corners_y, 0.0)
Xr = _poly_eval(calibration.fwd_a_r, corners_x, corners_y, 0.0)
Yr = _poly_eval(calibration.fwd_b_r, corners_x, corners_y, 0.0)
x_min = max(Xl.min(), Xr.min()) + margin_mm
x_max = min(Xl.max(), Xr.max()) - margin_mm
y_min = max(Yl.min(), Yr.min()) + margin_mm
y_max = min(Yl.max(), Yr.max()) - margin_mm
if x_min >= x_max or y_min >= y_max:
raise ValueError(
"Left and right cameras have no overlapping field of view "
f"(X: [{x_min:.1f}, {x_max:.1f}] Y: [{y_min:.1f}, {y_max:.1f}]). "
"Check that both cameras cover the same measurement plane.")
if scale_px_per_mm is None:
x_extent = ((Xl.max() - Xl.min()) + (Xr.max() - Xr.min())) / 2
scale_px_per_mm = w / x_extent
nx = int(round((x_max - x_min) * scale_px_per_mm))
ny = int(round((y_max - y_min) * scale_px_per_mm))
# world_y runs top-to-bottom: y_max at row 0 (physically high) down to
# y_min at the last row, so the backprojected image is right-side-up and
# px_size_y is negative (downward pixel → lower world Y).
world_x = np.linspace(x_min, x_max, nx)
world_y = np.linspace(y_max, y_min, ny)
return world_x, world_y
[docs]
def backproject_image(img, inv_a, inv_b, world_x, world_y, method='linear',
Z=0.0):
"""Back-project a raw camera image onto a regular world-coordinate grid.
For each point on the output world grid, the world->pixel mapping
(inv_a, inv_b) maps it to a source pixel location, which is then
interpolated from the raw image.
Parameters
----------
img : ndarray, shape (ny_img, nx_img)
Raw camera image (float or int).
inv_a, inv_b : ndarray
World->pixel mapping coefficients (pass calibration.map_a_l /
map_b_l etc.). The full 9-term Soloff map is evaluated at *Z*.
world_x : ndarray, shape (nx,)
X world coordinates of the output grid columns (mm), 1-D, monotone.
world_y : ndarray, shape (ny,)
Y world coordinates of the output grid rows (mm), 1-D.
method : str
Interpolation method passed to RegularGridInterpolator.
Z : float or ndarray, optional
Out-of-plane plane (mm) at which to evaluate the Soloff map. ``0.0``
(default) dewarps at the calibration plane (the in-plane Z=0 part of
the map, as before). A scalar shifts the whole grid to one Z; an array
must broadcast to the world grid ``(len(world_y), len(world_x))`` and
gives a per-point Δz — pass a Δz(X,Y) field to re-dewarp on the laser
sheet (disparity self-calibration iteration). Because the 9-term map
already carries the Z dependence, no extra geometry is needed.
Returns
-------
img_bp : ndarray, shape (ny, nx)
Back-projected image in world coordinates.
"""
from scipy.interpolate import RegularGridInterpolator
ny_img, nx_img = img.shape
ix = np.arange(nx_img, dtype=float)
iy = np.arange(ny_img, dtype=float)
interp = RegularGridInterpolator(
(ix, iy), img.T.astype(float),
method=method, bounds_error=False, fill_value=0.0,
)
Xg, Yg = np.meshgrid(world_x, world_y)
x_src = _poly_eval(inv_a, Xg, Yg, Z)
y_src = _poly_eval(inv_b, Xg, Yg, Z)
pts = np.stack([x_src.ravel(), y_src.ravel()], axis=1)
return interp(pts).reshape(Xg.shape)
[docs]
def list_images(src_dir, exts=IMAGE_EXTS):
"""Return the sorted image filenames in *src_dir* (any of *exts*)."""
return sorted(f for f in os.listdir(src_dir)
if f.lower().endswith(tuple(e.lower() for e in exts)))
[docs]
def download_pivchallenge_E(dest_dir, url=PIVCHALLENGE_E_URL,
force=False, keep_zip=False):
"""Download the 4th PIV Challenge (2014) Case E stereo dataset into
*dest_dir*, sorted into ``calibration``, ``camera1`` and ``camera2``
subfolders that hold *only* the image files.
Case E is a time-resolved stereo-PIV vortex-ring experiment recorded by
two cameras (1 and 3). The archive identifies every image by filename:
* calibration target ``E_camera_{1,3}_z_{1..7}.tif`` (7 Z-planes,
ΔZ = 1 mm) -> ``calibration`` (filenames keep the camera/plane, so both
cameras live in the same folder without ambiguity)
* particle frames ``E_camera_{1,3}_frame_NNNNN.tif`` -> ``camera1``
(the lower camera index) and ``camera2`` (the higher one)
Classification is by filename (``_z_`` vs ``_frame_`` and the ``camera_N``
index), so the layout inside the zip is irrelevant. Everything that is not
an image (the instructions PDF, readmes, ...) is dropped. This lets the
stereo tests run against a real dataset without committing ~200 MB of
images to the repository.
Parameters
----------
dest_dir : str
Target folder; ``calibration``/``camera1``/``camera2`` are created
inside it.
url : str, optional
Source archive (defaults to the official pivchallenge.org URL).
force : bool, optional
Re-download and re-extract even if the subfolders already hold images.
keep_zip : bool, optional
Keep the downloaded ``.zip`` (in *dest_dir*) instead of deleting it.
Returns
-------
dict
``{'calibration': path, 'camera1': path, 'camera2': path}`` with the
number of images sorted into each printed to stdout.
"""
sub = {name: os.path.join(dest_dir, name)
for name in ('calibration', 'camera1', 'camera2')}
for path in sub.values():
os.makedirs(path, exist_ok=True)
if not force and all(list_images(path) for path in sub.values()):
print(f"Case E already present in {dest_dir} (use force=True to "
"re-download).")
return sub
zip_path = os.path.join(dest_dir, os.path.basename(url))
if force or not os.path.isfile(zip_path):
print(f"Downloading Case E (~197 MB) from {url} ...")
# The host rejects the default urllib agent; spoof a browser one.
req = urllib.request.Request(url, headers={'User-Agent': 'Mozilla/5.0'})
tmp = zip_path + '.part'
with urllib.request.urlopen(req) as resp, open(tmp, 'wb') as out:
total = int(resp.headers.get('Content-Length', 0))
done = 0
while True:
chunk = resp.read(1 << 20)
if not chunk:
break
out.write(chunk)
done += len(chunk)
if total:
print(f"\r {done / 1e6:6.1f} / {total / 1e6:.1f} MB "
f"({100 * done / total:4.1f} %)", end='', flush=True)
print()
os.replace(tmp, zip_path)
raw_dir = os.path.join(dest_dir, '_case_E_raw')
if os.path.isdir(raw_dir):
shutil.rmtree(raw_dir)
print("Extracting ...")
with zipfile.ZipFile(zip_path) as zf:
zf.extractall(raw_dir)
cam_re = re.compile(r'camera_(\d+)', re.IGNORECASE)
# Map the two distinct camera indices (sorted) onto camera1 / camera2.
cameras = set()
for _, _, files in os.walk(raw_dir):
for fname in files:
m = cam_re.search(fname)
if m and fname.lower().endswith(IMAGE_EXTS):
cameras.add(int(m.group(1)))
cam_order = sorted(cameras)
cam_to_sub = {}
if len(cam_order) >= 1:
cam_to_sub[cam_order[0]] = sub['camera1']
if len(cam_order) >= 2:
cam_to_sub[cam_order[1]] = sub['camera2']
counts = {name: 0 for name in sub}
for root, _, files in os.walk(raw_dir):
for fname in files:
if not fname.lower().endswith(IMAGE_EXTS):
continue
src = os.path.join(root, fname)
if re.search(r'_z_\d+', fname, re.IGNORECASE):
dst_dir = sub['calibration']
key = 'calibration'
elif re.search(r'frame', fname, re.IGNORECASE):
m = cam_re.search(fname)
dst_dir = cam_to_sub.get(int(m.group(1))) if m else None
if dst_dir is None:
continue
key = os.path.basename(dst_dir)
else:
continue
shutil.move(src, os.path.join(dst_dir, fname))
counts[key] += 1
shutil.rmtree(raw_dir)
if not keep_zip:
os.remove(zip_path)
print(f"Case E ready in {dest_dir}: "
f"calibration={counts['calibration']}, "
f"camera1={counts['camera1']}, camera2={counts['camera2']} images.")
return sub
[docs]
def backproject_folder(src_dir, dst_dir, inv_a, inv_b, world_x, world_y,
exts=IMAGE_EXTS, files=None, Z=0.0):
"""Back-project images in src_dir onto the world-coordinate grid
and save results to dst_dir, preserving filenames and bit depth.
Parameters
----------
src_dir : str
Folder containing raw camera images (any format in *exts*).
dst_dir : str
Output folder (created if it does not exist).
inv_a, inv_b : ndarray
World->pixel mapping coefficients (pass calibration.map_a_l / map_b_l).
world_x, world_y : ndarray
World grid from make_world_grid().
exts : tuple of str
Image extensions to process (default :data:`IMAGE_EXTS`).
files : list of str or None
Explicit subset of filenames (relative to *src_dir*) to process. When
None (default) every image in the folder is dewarped; pass a slice of
``list_images(src_dir)`` to dewarp only selected frames.
Z : float or ndarray, optional
Out-of-plane plane (mm) for the dewarp, forwarded to
``backproject_image``. ``0.0`` (default) dewarps at the calibration
plane; pass a Δz(X,Y) field on the world grid to re-dewarp on the
laser sheet for disparity self-calibration iteration.
"""
os.makedirs(dst_dir, exist_ok=True)
for fname in (list_images(src_dir, exts) if files is None else files):
raw = cv2.imread(os.path.join(src_dir, fname), cv2.IMREAD_UNCHANGED)
img = rgb2gray(raw)
bp = backproject_image(img, inv_a, inv_b, world_x, world_y, Z=Z)
bp = np.clip(bp, 0, np.iinfo(raw.dtype).max).astype(raw.dtype)
cv2.imwrite(os.path.join(dst_dir, fname), bp)
[docs]
def compute_disparity(paths_cam1, paths_cam2,
window_size=128, overlap=0.75,
peak_search_ratio=0.6, peak_ratio=0.83,
show_corr_windows=False):
"""Compute the disparity field between two cameras from backprojected images.
Uses ensemble correlation: the normalised cross-correlation maps from all
image pairs are averaged before searching for the peak. This reinforces
the true disparity signal coherently while random noise averages toward
zero, giving a much sharper and more reliable peak than averaging
individual per-pair displacements.
Parameters
----------
paths_cam1 : list of str or ndarray
Back-projected world images for camera 1, given either as file paths
(any format ``cv2`` reads) or as pre-loaded arrays.
paths_cam2 : list of str or ndarray
Corresponding images for camera 2 (same order, same-instant pairs).
window_size : int
Interrogation window size in pixels (square). Default 128.
overlap : float
Fractional overlap between adjacent windows, in [0, 1). Default 0.75.
peak_search_ratio : float
Fraction of window_size used as the peak-search radius in the
correlation plane. Limits the maximum detectable disparity to
window_size * peak_search_ratio / 2 pixels. Default 0.6.
peak_ratio : float
Maximum allowed second-to-first correlation peak ratio (same
convention as DPIV.find_peaks). Windows whose secondary peak exceeds
this fraction of the primary are ambiguous and marked invalid (NaN).
Lower is stricter; >= 1 disables rejection. Default 0.83 (= 1/1.2).
show_corr_windows : bool
If True, display 4 sample mean correlation windows (spread across the
grid) with the detected peak marked, so you can inspect peak quality.
Default False.
Returns
-------
x_grid : ndarray, shape (ny_win, nx_win)
X pixel coordinate (column) of each window centre in the world image.
y_grid : ndarray, shape (ny_win, nx_win)
Y pixel coordinate (row) of each window centre in the world image.
disp_x : ndarray, shape (ny_win, nx_win)
X-disparity (cam2 minus cam1) in pixels of the world image, derived
from the ensemble-averaged correlation peak.
disp_y : ndarray, shape (ny_win, nx_win)
Y-disparity in pixels.
valid : ndarray of bool, shape (ny_win, nx_win)
True where the ensemble peak is unambiguous (peak_ratio not exceeded).
"""
if len(paths_cam1) != len(paths_cam2):
raise ValueError("paths_cam1 and paths_cam2 must have the same length.")
if not (0.0 <= overlap < 1.0):
raise ValueError("overlap must be in [0, 1).")
step = int(round(window_size * (1.0 - overlap)))
step = max(step, 1)
ws = window_size
half_ws = ws // 2
# --- determine grid from the first image pair ---
img1_ref = rgb2gray(paths_cam1[0])
ny_img, nx_img = img1_ref.shape
# window centres: first centre at half_ws, last where box fits inside image
cx = np.arange(half_ws, nx_img - half_ws + 1, step)
cy = np.arange(half_ws, ny_img - half_ws + 1, step)
x_grid, y_grid = np.meshgrid(cx, cy)
ny_win, nx_win = x_grid.shape
search_r = int(round(ws * peak_search_ratio / 2))
mid = ws // 2 # centre of correlation plane after fftshift
sz = 2 * search_r + 1
# Accumulate normalised cross-correlations across all pairs (ensemble).
# Each window gets its own (sz x sz) accumulator; dividing at the end
# gives the mean correlation whose peak is then found once.
corr_sum = np.zeros((ny_win, nx_win, sz, sz))
r_lo = mid - search_r
r_hi = mid + search_r + 1
c_lo = mid - search_r
c_hi = mid + search_r + 1
for path1, path2 in zip(paths_cam1, paths_cam2):
img1 = rgb2gray(path1)
img2 = rgb2gray(path2)
for jj in range(ny_win):
for ii in range(nx_win):
r0 = int(y_grid[jj, ii]) - half_ws
c0 = int(x_grid[jj, ii]) - half_ws
sub1 = img1[r0:r0 + ws, c0:c0 + ws]
sub2 = img2[r0:r0 + ws, c0:c0 + ws]
sub1 = sub1 - sub1.mean()
sub2 = sub2 - sub2.mean()
sig1 = max(1e-6, np.sqrt((sub1**2).sum()))
sig2 = max(1e-6, np.sqrt((sub2**2).sum()))
corr = np.fft.fftshift(np.abs(np.fft.ifft2(
np.multiply(np.conj(np.fft.fft2(sub1)),
np.fft.fft2(sub2))))) / (sig1 * sig2)
corr_sum[jj, ii] += corr[r_lo:r_hi, c_lo:c_hi]
# Mean ensemble correlation — find peak once per window
corr_mean = corr_sum / len(paths_cam1)
# Ambiguous windows stay NaN/invalid (no median_filter fallback here, unlike
# the main PIV pipeline) -> find_peaks is called with return_valid=True.
disp_x = np.full((ny_win, nx_win), np.nan)
disp_y = np.full((ny_win, nx_win), np.nan)
valid = np.zeros((ny_win, nx_win), dtype=bool)
for jj in range(ny_win):
for ii in range(nx_win):
cs = corr_mean[jj, ii]
eps_x, eps_y, mc, mr, ok = DPIV.find_peaks(
cs, sz, sz, westerweel=0,
peak_ratio=peak_ratio, return_valid=True)
if not ok:
continue
disp_x[jj, ii] = (mc + eps_x) - search_r
disp_y[jj, ii] = (mr + eps_y) - search_r
valid[jj, ii] = True
if show_corr_windows:
import matplotlib.pyplot as plt
# Pick 4 windows at quarter-points (25 % and 75 % along each axis)
# so they are equally spaced from both centre and border
rows_idx = [int(round(ny_win * q)) for q in (0.25, 0.75)]
cols_idx = [int(round(nx_win * q)) for q in (0.25, 0.75)]
# clamp to valid range
rows_idx = [max(0, min(ny_win - 1, r)) for r in rows_idx]
cols_idx = [max(0, min(nx_win - 1, c)) for c in cols_idx]
samples = [(jj, ii) for jj in rows_idx for ii in cols_idx]
px_axis = np.arange(-search_r, search_r + 1, dtype=float)
PX, PY = np.meshgrid(px_axis, px_axis) # PX = dx axis, PY = dy axis
fig, axes = plt.subplots(2, 2, figsize=(10, 8))
axes = axes.ravel()
for ax, (jj, ii) in zip(axes, samples):
cs = corr_mean[jj, ii]
cf = ax.contourf(PX, PY, cs, levels=30, cmap='plasma')
fig.colorbar(cf, ax=ax)
if valid[jj, ii]:
ax.plot(disp_x[jj, ii], disp_y[jj, ii],
'+', color='white', ms=12, mew=2)
ax.set_title(f'window ({ii}, {jj}) '
f'dx={disp_x[jj, ii]:.2f} dy={disp_y[jj, ii]:.2f}',
fontsize=9)
ax.set_xlabel('dx (px)', fontsize=8)
ax.set_ylabel('dy (px)', fontsize=8)
ax.set_aspect('equal')
fig.suptitle('Ensemble mean correlation — 4 sample windows', fontsize=12)
plt.tight_layout()
plt.show()
return x_grid, y_grid, disp_x, disp_y, valid
[docs]
def build_disparity_maps(x_d, y_d, disp_x, disp_y, valid, height, width):
"""Interpolate the sparse disparity field to full world-image resolution.
The sparse grid returned by ``compute_disparity`` has one vector per
interrogation window. This function lifts it to a dense (H x W) map
so that ``apply_disparity_correction`` can warp every pixel cheaply
with a single ``cv2.remap`` call. Call this once before the PIV loop
and reuse the result for every image.
Parameters
----------
x_d, y_d : ndarray, shape (ny_win, nx_win)
Pixel coordinates of the sparse disparity window centres.
disp_x, disp_y : ndarray, shape (ny_win, nx_win)
Disparity in world-image pixels (cam2 minus cam1).
valid : ndarray of bool, shape (ny_win, nx_win)
Mask of reliable vectors. Invalid entries are filled by
nearest-neighbour extrapolation before interpolating.
height, width : int
Size of the world image in pixels.
Returns
-------
dx_full, dy_full : ndarray, shape (H, W), dtype float32
Dense disparity maps ready for ``apply_disparity_correction``.
"""
from scipy.interpolate import griddata
pts_valid = np.column_stack((x_d[valid].ravel(), y_d[valid].ravel()))
dx_valid = disp_x[valid].ravel()
dy_valid = disp_y[valid].ravel()
gx = np.arange(width, dtype=np.float32)
gy = np.arange(height, dtype=np.float32)
Xg, Yg = np.meshgrid(gx, gy)
pts_full = np.column_stack((Xg.ravel(), Yg.ravel()))
dx_full = griddata(pts_valid, dx_valid, pts_full,
method='cubic', fill_value=np.nan).reshape(height, width)
dy_full = griddata(pts_valid, dy_valid, pts_full,
method='cubic', fill_value=np.nan).reshape(height, width)
# Fill border NaN (outside convex hull) with nearest-neighbour values
nan_mask = np.isnan(dx_full)
if nan_mask.any():
dx_nn = griddata(pts_valid, dx_valid, pts_full,
method='nearest').reshape(height, width)
dy_nn = griddata(pts_valid, dy_valid, pts_full,
method='nearest').reshape(height, width)
dx_full[nan_mask] = dx_nn[nan_mask]
dy_full[nan_mask] = dy_nn[nan_mask]
return dx_full.astype(np.float32), dy_full.astype(np.float32)
[docs]
def disparity_to_delta_z(calibration, x_grid, y_grid, disp_x, disp_y, valid,
world_x, world_y, img_shape, piv_shape,
surface_order=None):
"""Convert a measured stereo disparity field to an out-of-plane Δz map.
The disparity is the residual misregistration left when a quiescent,
same-instant pair is dewarped to the common world grid at Z=0 and
cross-correlated (``backproject_*`` + ``compute_disparity``). It is
non-zero because the particles live on the laser sheet, at an offset
Δz(X,Y) from the calibration plane (Wieneke 2005). The returned Δz is
sampled on the Soloff reconstruction grid (``_stereo_recon_grid``), ready
to pass to ``soloff_velocity(..., delta_z=...)``.
``compute_disparity`` returns the disparity in *world-image pixels* on
a *world-image-pixel* grid, whereas the map Jacobian's parallax slope
∂(x,y)/∂Z is in *raw-camera pixels per mm*. These must be reconciled:
1. Window-centre pixel coords (x_grid, y_grid) and the disparity are
converted to world millimetres via the world-image scale (the spacing
of ``world_x``/``world_y`` used for the back-projection).
2. A particle at Z=Δz, dewarped assuming Z=0, lands at world offset
J⁻¹·(∂map/∂Z)·Δz in that camera's dewarped image (J = in-plane map
Jacobian; J⁻¹ turns the raw-pixel parallax into the world-mm shift the
dewarp actually shows). The cam2-minus-cam1 disparity is therefore
d_world = G · Δz, G = J_r⁻¹·∂map_r/∂Z − J_l⁻¹·∂map_l/∂Z
(a dimensionless 2-vector). Δz is recovered by least squares
(projecting the 2-D disparity onto G): Δz = (d_world·G)/(G·G).
Sign convention: cam1 must be the LEFT camera (``map_*_l``) and cam2 the
RIGHT (``map_*_r``), matching ``compute_disparity``'s cam2-minus-cam1
definition. **Qualitatively validated** -- iterating the self-cal loop
(``stereo_tutorial.py``) for 2 and 3 passes improves all metrics (the residual
disparity converges toward center), which confirms the Δz sign and a self-
consistent magnitude. Still to do: check the *absolute* Δz/w against an
independent ground truth (e.g. Willert's pivbook vortex ring or a synthetic
known-``w`` field).
Parameters
----------
calibration : Parameters.Stereo_Calibration
Must carry the dual-plane Soloff maps ``map_a_l/map_b_l/map_a_r/map_b_r``.
x_grid, y_grid, disp_x, disp_y, valid : ndarray
Outputs of ``compute_disparity`` (window-centre pixel coords, disparity
in world-image pixels, and the reliability mask).
world_x, world_y : ndarray, 1-D
World axes (mm) of the back-projected images the disparity was measured
on (the ``make_world_grid`` arrays used for back-projection). Their
spacing sets the world-image pixel scale.
img_shape, piv_shape : tuple
Raw image size and PIV grid shape, forwarded to ``_stereo_recon_grid``
so the returned Δz lands on the same grid ``soloff_velocity`` uses.
surface_order : None, 1 or 2, optional
Optional low-order surface fit of the scattered per-point Δz before it is
sampled onto the reconstruction grid -- a regularizer that denoises the
Δz field (a laser sheet's Δz is physically low-order: a flat/tilted/gently
bowed plane). The raw per-point estimate ``Δz=(d·G)/(G·G)`` is
noise-limited (its error tracks the disparity noise ~1:1); fitting a
surface averages that noise out.
* ``None`` (default) -- raw per-point + ``griddata`` resample (the
original behaviour; existing runs are unchanged).
* ``1`` -- fit a plane ``[1, X, Y]`` (≈ Wieneke's exact triangulation).
* ``2`` -- fit a quadratic ``[1, X, Y, X², XY, Y²]``; keeps real sheet
curvature *and* rejects noise (best for a bowed sheet, where a plane
is biased). Going higher would re-admit noise -- don't.
A synthetic geometric study (iid disparity noise) found order 1/2 beat
the raw per-point Δz ~30× on flat/tilted sheets, and on a curved sheet
order 2 beats both the raw estimate (~20×) and a plane (which is frozen
at the curvature bias). NB: this improves *noise behaviour only* -- it
is not a substitute for validating the *absolute* Δz against ground truth.
Returns
-------
delta_z : ndarray, shape == reconstruction grid (= piv_shape)
Out-of-plane offset (mm) per grid point, pluggable as
``soloff_velocity(..., delta_z=delta_z)``.
"""
from scipy.interpolate import griddata
maps = (calibration.map_a_l, calibration.map_b_l,
calibration.map_a_r, calibration.map_b_r)
if any(c is None or np.any(np.isnan(np.asarray(c, dtype=float)))
for c in maps):
raise ValueError("disparity_to_delta_z requires dual-plane Soloff maps "
"(map_a_l/map_b_l/map_a_r/map_b_r missing or NaN).")
map_a_l, map_b_l, map_a_r, map_b_r = maps
# World-image pixel -> world mm. world_y runs top-to-bottom (decreasing),
# so px_size_y is negative; both disp and centre coords use the same scale.
W, H = len(world_x), len(world_y)
px_x = (world_x[-1] - world_x[0]) / (W - 1)
px_y = (world_y[-1] - world_y[0]) / (H - 1)
m = valid.ravel()
xc = x_grid.ravel()[m]
yc = y_grid.ravel()[m]
Xw = world_x[0] + xc * px_x # window centres in world mm
Yw = world_y[0] + yc * px_y
d_mm_x = disp_x.ravel()[m] * px_x # disparity in world mm (cam2 - cam1)
d_mm_y = disp_y.ravel()[m] * px_y
def _Jinv_dZ(map_a, map_b):
"""J^{-1}·(∂map/∂Z) at Z=0 for one camera: the world-mm apparent shift
per mm of Z that survives the Z=0 dewarp. Returns (gx, gy) arrays."""
ax, ay, az = _poly_jacobian(map_a, Xw, Yw) # ∂x/∂X, ∂x/∂Y, ∂x/∂Z
bx, by, bz = _poly_jacobian(map_b, Xw, Yw) # ∂y/∂X, ∂y/∂Y, ∂y/∂Z
det = ax * by - ay * bx # det(J), per point
gx = (by * az - ay * bz) / det # (J^{-1}·[az,bz])_x
gy = (-bx * az + ax * bz) / det # (J^{-1}·[az,bz])_y
return gx, gy
gxl, gyl = _Jinv_dZ(map_a_l, map_b_l)
gxr, gyr = _Jinv_dZ(map_a_r, map_b_r)
Gx = gxr - gxl # G = J_r^{-1}dZ_r - J_l^{-1}dZ_l
Gy = gyr - gyl
denom = Gx * Gx + Gy * Gy
dz_pts = (d_mm_x * Gx + d_mm_y * Gy) / denom
Xg, Yg = _stereo_recon_grid(calibration, img_shape, piv_shape)
if surface_order is None:
# Resample the raw scattered Δz onto the Soloff reconstruction grid.
pts = np.column_stack((Xw, Yw))
dz = griddata(pts, dz_pts, (Xg, Yg), method='linear')
nan = np.isnan(dz)
if nan.any():
dz[nan] = griddata(pts, dz_pts, (Xg, Yg), method='nearest')[nan]
return dz
if surface_order not in (1, 2):
raise ValueError(
f"surface_order must be None, 1 or 2 (got {surface_order!r}).")
# Low-order surface fit of the scattered per-point Δz, evaluated analytically
# on the grid (no griddata). Normalise coords by the grid extent so the fit
# is well conditioned and the scattered points and grid share one frame.
x0, xr = Xg.min(), max(Xg.ptp(), 1e-12)
y0, yr = Yg.min(), max(Yg.ptp(), 1e-12)
def _design(xn, yn):
"""Low-order surface-fit design matrix on normalized coords:
[1, X, Y] for order 1, plus [X², XY, Y²] for order 2."""
cols = [np.ones_like(xn), xn, yn]
if surface_order == 2:
cols += [xn * xn, xn * yn, yn * yn]
return np.stack([c.ravel() for c in cols], axis=1)
A = _design((Xw - x0) / xr, (Yw - y0) / yr)
coef, *_ = np.linalg.lstsq(A, dz_pts, rcond=None)
dz = (_design((Xg - x0) / xr, (Yg - y0) / yr) @ coef).reshape(Xg.shape)
return dz
[docs]
def delta_z_on_world_grid(delta_z, world_x, world_y):
"""Resample a recon-grid Δz field onto the world-image dewarp grid.
``disparity_to_delta_z`` returns Δz on the Soloff *reconstruction* grid
(``piv_shape``), but re-dewarping for disparity self-calibration iteration
needs Δz on the *world-image* grid (``world_x`` x ``world_y``, full image
resolution), so it can be passed straight as ``Z`` to ``backproject_*``.
The reconstruction grid is just ``world_x``/``world_y`` resampled to PIV
resolution (see ``_stereo_recon_grid``), so the two grids share the same
world rectangle and the recon-grid axes are recovered from the world
endpoints and ``delta_z``'s shape -- no calibration or image size needed.
``delta_z`` may also be ``None`` (-> ``0.0``: dewarp at the calibration
plane) or a scalar (uniform offset); both pass straight through, so the
caller can hand whatever Δz estimate it has to ``backproject_*(..., Z=...)``
with no branch.
Parameters
----------
delta_z : None, float, or ndarray
Δz (mm): ``None`` / scalar / recon-grid field (shape = ``piv_shape``).
world_x, world_y : ndarray, 1-D
World axes (mm) of the dewarp grid (the ``make_world_grid`` arrays).
Returns
-------
Z_world : float or ndarray, shape (len(world_y), len(world_x))
Δz on the world grid, ready as ``backproject_image(..., Z=Z_world)``.
"""
if delta_z is None:
return 0.0
dz = np.asarray(delta_z, dtype=float)
if dz.ndim == 0:
return float(dz)
ny, nx = dz.shape
rx = np.linspace(world_x[0], world_x[-1], nx)
ry = np.linspace(world_y[0], world_y[-1], ny)
Xg, Yg = np.meshgrid(rx, ry)
src = np.column_stack((Xg.ravel(), Yg.ravel()))
Xw, Yw = np.meshgrid(world_x, world_y)
out = griddata(src, dz.ravel(), (Xw, Yw), method='linear')
nan = np.isnan(out)
if nan.any():
out[nan] = griddata(src, dz.ravel(), (Xw, Yw), method='nearest')[nan]
return out
[docs]
def disparity_convergence_report(disp_x, disp_y, valid, dz_inc, prev_rms=None,
tag=""):
"""Print disparity self-calibration metrics; return disp_rms for the next call.
disp_rms alone plateaus at the PIV noise floor, so also show: the Δz step
RMS (-> 0 when converged), the bias (coherent part Δz removes -> 0), the
scatter (~irreducible floor) and the valid-window count. A disp_rms that
GROWS vs ``prev_rms`` means a wrong Δz sign for this geometry.
"""
m = np.asarray(valid, dtype=bool)
dx, dy = disp_x[m], disp_y[m]
disp_rms = np.sqrt(np.nanmean(dx**2 + dy**2))
bias = np.hypot(np.nanmean(dx), np.nanmean(dy))
scatter = np.sqrt(np.nanvar(dx) + np.nanvar(dy))
step = np.sqrt(np.nanmean(np.asarray(dz_inc, float)**2))
trend = ""
if prev_rms is not None:
trend = " GREW - check Δz sign!" if disp_rms > prev_rms * 1.05 else " (down)"
prefix = f"{tag} " if tag else " "
print(f"{prefix}Δz step={step:.4f} mm | disp_rms={disp_rms:.3f} px{trend} | "
f"bias={bias:.3f} | noise={scatter:.3f} | valid={int(m.sum())}")
return disp_rms
[docs]
def calibration_GUI(show_diagnostics=False):
"""Launch the stereo calibration GUI (defined in stereo_gui.py).
Imported lazily so this math module has no hard Tkinter/PIL dependency.
Parameters
----------
show_diagnostics : bool
Developer flag, forwarded to the GUI: when True, ``sort_points`` pops
its diagnostic plot (the detected dots annotated with their assigned
lattice (row, col) labels). Off by default so the workflow is not
interrupted.
"""
from dpivsoft.stereo_gui import calibration_GUI as _launch
_launch(show_diagnostics=show_diagnostics)
if __name__ == "__main__":
calibration_GUI()