Stereo PIV: calibration & 3-component reconstruction¶

This page explains the stereo-PIV calibration and 3-component (\(u, v, w\)) reconstruction implemented in dpivsoft/stereo_DPIV.py (calibration GUI in dpivsoft/stereo_gui.py): the two polynomial calibration models and their fitting, the local-Jacobian Soloff reconstruction, and the disparity-based \(\Delta z\) self-calibration.

1. The core idea¶

Two cameras view the light sheet from different angles, so a particle moving out-of-plane (\(W\)) appears in each camera as an apparent in-plane image shift whose size and direction depend on that camera’s viewing geometry. The Soloff method inverts this: at each grid point it asks what world motion \((U, V, W)\) best explains the in-plane image shifts both cameras actually measured.

Answering that requires, per camera, the local derivative of the world→pixel mapping — its Jacobian — which the calibration polynomial provides analytically at every point in the field. Reconstruction runs directly on the raw camera images: no global viewing angle, no dewarping in the reconstruction path. Because the viewing geometry is recovered locally, the method stays accurate across the whole field of view.

2. Coordinate conventions¶

World coordinates \((X, Y, Z)\) in mm. \(Z = 0\) is the calibration / light-sheet plane; \(Z\) is the out-of-plane axis and \(W\) the out-of-plane velocity component.
Image coordinates \((x, y)\) per camera, in pixels. Subscripts \(l\) / \(r\) denote the left / right cameras.
Camera configuration: Front-Front (both cameras on the same side of the sheet) or Front-Back (opposite sides). In Front-Back the left camera images the sheet from behind, so get_real_* flips the left camera’s world \(X\) — the geometry is baked into the fitted maps and reconstruction needs no extra sign handling. The configuration is stored with the calibration.
Sign of \(w\): positive \(w\) points toward the right camera, negative toward the left.

3. The two calibration models¶

Each camera’s map is a polynomial in \((p, q, z)\), where the meaning of \(p, q\) depends on the fit direction (§4). The model is selected by the mode argument of stereo_calibration() (and in the GUI’s “Soloff Calculation” popup); consumers auto-detect the model of a stored calibration from the coefficient-vector length (9 vs 19), so a saved calibration always carries its own model.

3.1 Reduced — 9 terms (`mode='reduced'`, default)¶

Quadratic in-plane at order \(z^0\), linear in-plane at order \(z^1\):

\[ F(p, q, z) \;=\; c_0\,p^2 + c_1\,q^2 + c_2\,pq + c_3\,p + c_4\,q + c_5 \;+\; \left(c_6 + c_7\,p + c_8\,q\right) z \]

The last bracket is the parallax block — the only depth-dependent part. Constraining it requires dots at ≥ 2 distinct \(Z\) planes. Design matrix: _poly9_design(); fit: _fit_poly9().

3.2 Full — 19 terms (`mode='full'`), the classic Soloff model¶

Cubic in-plane at \(z^0\), quadratic at \(z^1\), linear at \(z^2\) (Soloff et al. 1997):

\[\begin{split} \begin{aligned} F(p, q, z) \;=\;& c_0\,p^3 + c_1\,p^2q + c_2\,pq^2 + c_3\,q^3 + c_4\,p^2 + c_5\,pq + c_6\,q^2 + c_7\,p + c_8\,q + c_9 && (z^0)\\[2pt] &+ \left(c_{10}\,p^2 + c_{11}\,pq + c_{12}\,q^2 + c_{13}\,p + c_{14}\,q + c_{15}\right) z && (z^1)\\[2pt] &+ \left(c_{16}\,p + c_{17}\,q + c_{18}\right) z^2 && (z^2) \end{aligned} \end{split}\]

The cubic in-plane terms \(p^3, p^2q, pq^2, q^3\) model lens distortion (radial distortion is leading-order \(r^3\), so nothing below cubic can represent it). The asymmetric orders are deliberate: many in-plane dots determine the cubic terms well, while few depth planes limit the model to quadratic in \(z\).

Requirements: ≥ 3 distinct \(Z\) planes (the highest \(z\)-power) and enough in-plane dots per plane for the cubic terms. Design matrix: _poly19_design(); fit: _fit_poly19(). The GUI refuses to run the full model with fewer than 3 distinct \(Z\) planes per camera.

4. Fitting the calibration (`stereo_calibration`)¶

For each camera, two maps are fitted by exact least squares (np.linalg.lstsq):

fwd_* — pixel → world: \((X, Y) = \mathrm{poly}(x_\mathrm{px}, y_\mathrm{px}, Z)\). Its \(Z = 0\) part builds the world-coordinate grid (make_world_grid()).
map_* — world → pixel (the Soloff map): \((x_\mathrm{px}, y_\mathrm{px}) = \mathrm{poly}(X, Y, Z)\). Used by soloff_velocity() (local Jacobian) and by back-projection.

So each camera has four coefficient vectors — fwd_a, fwd_b (world \(X, Y\)) and map_a, map_b (pixel \(x, y\)) — each a 9-vector (reduced) or 19-vector (full). A dual-plane calibration (≥ 2 distinct \(Z\)) is the minimum that constrains the parallax; if the plane requirement is not met, _fit_camera returns NaN maps and the result reports “Soloff maps: NOT available”.

4.1 Conditioning¶

The raw design columns span many orders of magnitude (\(p^2\) in pixels \(\sim 10^7\) vs the \(z\) column \(\sim 1\)): the fit is badly conditioned and the weakly-excited parallax coefficients — the ones stereo depends on — would drown in noise. Two equivalent-in-goal strategies:

Reduced (_fit_poly9): fit in normalised coordinates (each of \(p, q, z\) shifted to zero mean, unit standard deviation), then fold the coefficients back into the raw basis algebraically — the 9-term basis is closed under per-axis affine normalisation. Condition number drops from \(\sim 10^7\) to \(\sim 1\).
Full (_fit_poly19): column equilibration — scale each design column to unit RMS before the solve, divide the coefficients back afterwards. Leaves the least-squares solution unchanged, and is basis-agnostic, so it covers the cubic / \(z^2\) columns without a fold-back derivation.

Both return coefficients in the raw basis: downstream evaluators and the on-disk format depend only on the coefficient length, not on the conditioning.

5. The local mapping Jacobian¶

The analytic derivatives of the world→pixel map_* polynomials give how an image point moves per unit of world motion. Writing \(F\) for one image coordinate (\(x\) from map_a, \(y\) from map_b):

Reduced (9-term), at any \(Z\):

\[\begin{split} \begin{aligned} \frac{\partial F}{\partial X} &= 2c_0 X + c_2 Y + c_3 + c_7 Z\\[2pt] \frac{\partial F}{\partial Y} &= 2c_1 Y + c_2 X + c_4 + c_8 Z\\[2pt] \frac{\partial F}{\partial Z} &= c_6 + c_7 X + c_8 Y \qquad\text{(parallax slope; independent of } Z\text{)} \end{aligned} \end{split}\]

Full (19-term), at any \(Z\):

\[\begin{split} \begin{aligned} \frac{\partial F}{\partial X} &= 3c_0 X^2 + 2c_1 XY + c_2 Y^2 + 2c_4 X + c_5 Y + c_7 + 2c_{10} XZ + c_{11} YZ + c_{13} Z + c_{16} Z^2\\[2pt] \frac{\partial F}{\partial Y} &= c_1 X^2 + 2c_2 XY + 3c_3 Y^2 + c_5 X + 2c_6 Y + c_8 + c_{11} XZ + 2c_{12} YZ + c_{14} Z + c_{17} Z^2\\[2pt] \frac{\partial F}{\partial Z} &= c_{10} X^2 + c_{11} XY + c_{12} Y^2 + c_{13} X + c_{14} Y + c_{15} + 2c_{16} XZ + 2c_{17} YZ + 2c_{18} Z \end{aligned} \end{split}\]

In the full model the parallax slope \(\partial F/\partial Z\) is itself \(Z\)-dependent (the curvature the third plane buys), and the in-plane derivatives carry the cubic-distortion terms.

In both models \(\partial(x, y)/\partial Z\) — the parallax — is what makes \(w\) observable, and it varies across the field (e.g. the \(c_7 X + c_8 Y\) part): the effective viewing angle is local, not a single global number. The derivatives default to \(Z = 0\) (the calibration plane); pass a \(\Delta z\) field for disparity-corrected reconstruction (§10).

6. Soloff 3C reconstruction (`soloff_velocity`)¶

Reconstruction happens on a common world grid at PIV vector resolution (not full image resolution). At each grid point:

Map the grid point into each camera’s raw image (_poly_eval() on map_*) and interpolate that camera’s measured PIV displacement there (griddata).
Build each camera’s Jacobian (_poly_jacobian() on map_*).

For camera \(i\), the Jacobian links world displacement to that camera’s measured image displacement:

\[\begin{split} \begin{bmatrix} \Delta x_i \\ \Delta y_i \end{bmatrix} = \begin{bmatrix} \dfrac{\partial x}{\partial X} & \dfrac{\partial x}{\partial Y} & \dfrac{\partial x}{\partial Z}\\[6pt] \dfrac{\partial y}{\partial X} & \dfrac{\partial y}{\partial Y} & \dfrac{\partial y}{\partial Z} \end{bmatrix} \begin{bmatrix} \Delta X \\ \Delta Y \\ \Delta Z \end{bmatrix} \end{split}\]

Stacking both cameras gives 4 equations for 3 unknowns, solved by least squares at every grid point:

\[\begin{split} \mathbf{d} = \begin{bmatrix} \Delta x_l \\ \Delta y_l \\ \Delta x_r \\ \Delta y_r \end{bmatrix} = M \begin{bmatrix} \Delta X \\ \Delta Y \\ \Delta Z \end{bmatrix} \qquad\Longrightarrow\qquad (U, V, W) = \left(M^{\top} M\right)^{-1} M^{\top}\, \mathbf{d} \end{split}\]

The two cameras’ parallaxes have opposite sign, so their disagreement isolates \(W\) while their agreement gives the in-plane \(U, V\). Grid points where either camera’s interpolated displacement is invalid (outside its field of view) are masked to NaN.

7. Model dispatch¶

The polynomial maths lives in explicit per-model functions (_poly9_design / _poly9_eval / _poly9_jacobian / _fit_poly9, and the _poly19_* counterparts) behind two thin selectors, _poly_eval and _poly_jacobian, which dispatch on coefficient length via _is_full_poly() — length 19 → full, 9 → reduced, anything else raises (a malformed map is an error, never a silent fallback). Every consumer — soloff_velocity, backproject_image, disparity_to_delta_z, make_world_grid, the calibration diagnostics — goes through the selectors, so a loaded 9- or 19-length calibration automatically uses the correct model.

8. When does the full model help?¶

Even evaluated at \(Z = 0\), where the \(z^2\) and cubic terms drop out, the full model is a different fit: the polynomial columns are correlated, so the extra columns change the other coefficients. The improvement scales with how much real \(Z\)-curvature and lens distortion the optics have over the calibration depth:

Parallax (→ \(W\)): with symmetric \(Z\) planes the reduced model already recovers the parallax slope at the centre; with asymmetric spacing the reduced slope is a chord across the \(Z\) range rather than the tangent at the sheet — the full model fixes that. Off the calibration plane (\(Z = \Delta z \neq 0\), the disparity path) the \(z^2\) terms are directly active.
In-plane (→ \(U, V\), positions): the reduced model aliases the mean of the un-modelled curvature into its constant / in-plane terms; the full model removes that bias.
Lens distortion: only the cubic in-plane terms can represent it at all.

A lower fit RMS is not automatically a better measurement: with exactly 3 planes the quadratic-in-\(Z\) fit passes through the three plane means exactly (no redundancy in \(Z\)), so part of any RMS drop is fitted noise. Prefer 4+ planes and validate with leave-one-plane-out reprojection or a synthetic known-\(W\) test. For thin sheets, near-telecentric or low-distortion optics, the gain is negligible.

9. World grid & dewarp¶

make_world_grid maps the four image corners of both cameras to world coordinates at \(Z = 0\) (_poly_eval on fwd_*), intersects them to the cameras’ overlap region, and returns a world axis grid (mm) — at full image resolution for dewarping, or sampled to the PIV shape inside _stereo_recon_grid.
backproject_image() dewarps a raw image onto a world grid by evaluating a world→pixel map (_poly_eval) to find the source pixel for each world point, then interpolating. Both work for 9- and 19-length maps via the selectors.

10. Disparity / \(\Delta z\) self-calibration¶

An optional refinement: by default reconstruction runs at \(Z = 0\); passing a \(\Delta z\) field to soloff_velocity enables the correction.

10.0 The problem¶

Calibration declares the target plane to be \(Z = 0\) and fits the maps there. In the experiment the particles are lit by the laser sheet, which in general sits at an offset — possibly tilted or curved — \(Z = \Delta z(X, Y)\) from that plane. Reconstructing as if the particles were at \(Z = 0\) evaluates the map, parallax included, at the wrong depth; since parallax is the only source of depth information, this biases \(w\) first and \((u, v)\) second. Self-calibration (Wieneke 2005) measures \(\Delta z\) from the particle images themselves and re-evaluates the calibration at \(Z = \Delta z\).

10.1 Step 1 — measure the disparity (Willert 1997)¶

Take a simultaneous frame pair — both cameras at the same instant, so they photograph the physically same particles and flow contributes nothing. Dewarp each frame onto the common world grid as if every particle were at \(Z = 0\) (§9), then cross-correlate the two dewarped images:

Sheet on the calibration plane → the same particle lands on the same world spot in both frames → zero disparity.
Sheet at \(Z = \Delta z\) → the cameras’ different viewing angles drop the particle on different world spots; the correlation measures the residual gap — the disparity field \(\mathbf{d}(X, Y) = (d_x, d_y)\), defined right-camera minus left-camera.

compute_disparity() returns \(\mathbf{d}\) in world-image pixels; it is converted to world millimetres for Step 2.

10.2 Step 2 — forward model: disparity per unit \(\Delta z\)¶

Split each camera’s world→pixel Jacobian (§5) into the in-plane \(2 \times 2\) block and the parallax column:

\[\begin{split} \frac{\partial (x, y)}{\partial (X, Y, Z)} \;=\; \Big[\;\; J \;\;\Big|\;\; \partial M/\partial Z \;\;\Big], \qquad J = \begin{bmatrix} a_x & a_y \\ b_x & b_y \end{bmatrix}, \quad \partial M/\partial Z = (a_z, b_z) \end{split}\]

A particle at true depth \(\Delta z\), dewarped assuming \(Z = 0\), leaves in camera \(i\) a raw-pixel error \((\partial M_i/\partial Z)\,\Delta z\). The disparity, however, is measured on the dewarped world grid in millimetres, so the pixel error must be pulled back through the in-plane inverse: camera \(i\)’s apparent world-mm shift is

\[ \mathbf{s}_i \;=\; J_i^{-1}\, \frac{\partial M_i}{\partial Z}\, \Delta z \]

(computed by _Jinv_dZ). Without the \(J^{-1}\) factor a pixel quantity would be compared with a millimetre one. Subtracting right minus left, consistent with Step 1:

\[ \mathbf{d} \;=\; \mathbf{s}_r - \mathbf{s}_l \;=\; \left( J_r^{-1} \frac{\partial M_r}{\partial Z} \;-\; J_l^{-1} \frac{\partial M_l}{\partial Z} \right) \Delta z \;\equiv\; \mathbf{G}\, \Delta z \]

\(\mathbf{G}(X, Y)\) is the world-mm disparity produced per mm of \(\Delta z\) — the net parallax that survives a \(Z = 0\) dewarp. The two cameras’ parallaxes point in opposite directions, so the subtraction adds their magnitudes: \(\mathbf{G}\) is large where the viewing directions differ strongly and shrinks toward zero where the cameras share a line of sight — there \(\Delta z\) is geometrically unobservable. \(\mathbf{G}\) is recomputed at every grid point from _poly_jacobian (a field, not a global constant; for the 19-term map \(\partial M/\partial Z\) is itself \(Z\)-dependent, which the same code handles). The pixel→world fwd_* fit is not used in this path.

10.3 Step 3 — least-squares inversion¶

\(\mathbf{d} = \mathbf{G}\,\Delta z\) is two equations (\(x\) and \(y\)) for one unknown, and \(\mathbf{d}\) is noisy. The least-squares solution is the projection

\[ \Delta z \;=\; \frac{\mathbf{d} \cdot \mathbf{G}}{\mathbf{G} \cdot \mathbf{G}} \;=\; \frac{d_x G_x + d_y G_y}{G_x^2 + G_y^2} \]

which keeps the component of \(\mathbf{d}\) along \(\mathbf{G}\) and discards the orthogonal component — under the model \(\mathbf{d} = \mathbf{G}\,\Delta z\), that component can only be noise.

This is the classic stereo formula \(\Delta z = \Delta u / (\tan\alpha_1 + \tan\alpha_2)\) made local and 2-D: \(J_i^{-1}\,\partial M_i/\partial Z\) is the local \(\tan\alpha_i\) computed point-by-point instead of assumed constant, subtracting the opposite-side cameras is the sum of tangents, and the projection is the least-squares version of the single division.

Surface regularisation (surface_order). The raw per-point projection passes disparity noise straight through into \(\Delta z\). Optionally, disparity_to_delta_z() fits a low-order surface to the scattered per-point \(\Delta z\) — basis \([1, X, Y]\) for order 1, plus \([X^2, XY, Y^2]\) for order 2 — and evaluates the fit analytically on the reconstruction grid:

surface_order=None — raw per-point field (default);
surface_order=1 — a plane, equivalent to Wieneke’s planar self-calibration without adopting his pinhole model;
surface_order=2 — keeps sheet tilt and curvature while averaging out per-point noise; the best performer in synthetic comparisons.

10.4 Applying the correction¶

soloff_velocity(..., delta_z=...) evaluates the map and Jacobian at \(Z = \Delta z(X, Y)\) instead of \(Z = 0\) — no image re-warp, because the calibration polynomial already carries the \(Z\) dependence. Re-measuring the disparity after correction should shrink it; iterating drives it toward zero and doubles as a built-in quality check.

10.5 Assumptions¶

First order in \(\Delta z\): \(\mathbf{s}_i = J_i^{-1} (\partial M_i/ \partial Z)\, \Delta z\) linearises the map about \(Z = 0\) — valid for a small sheet-to-target offset; the 19-term map’s \(Z\)-dependence and the iteration recover higher-order effects.
Simultaneity: the frame pair must be same-instant so the disparity is pure geometry, not flow. \(\Delta z\) is treated as quasi-static (a one-time calibration step).
Noise model: the projection is optimal under isotropic Gaussian disparity noise and degrades gracefully where \(\mathbf{G}\) is weak.

11. File format (`.npz`)¶

Saved by the GUI:

fwd_a_l, fwd_b_l, fwd_a_r, fwd_b_r     # pixel→world coeffs (X, Y) per camera
map_a_l, map_b_l, map_a_r, map_b_r     # world→pixel coeffs (x, y) per camera
camera_config                          # 'Front-Front' | 'Front-Back'

Each coefficient vector has shape (9,) (reduced) or (19,) (full). The loader (in Classes.py) detects the model from the coefficient length, so both calibrations load and run through the same path.

12. Function reference¶

function	role
`stereo_calibration()`	fit both cameras, both directions; returns the 8 coefficient vectors and a results string
`_fit_poly9()` / `_fit_poly19()`	least-squares fit of one map (reduced / full)
`_poly9_design()` / `_poly19_design()`	design matrices (also used for the fit RMS)
`_poly9_eval()` / `_poly19_eval()` / `_poly_eval()`	evaluate a map (reduced / full / dispatch)
`_poly9_jacobian()` / `_poly19_jacobian()` / `_poly_jacobian()`	map derivatives (reduced / full / dispatch)
`_is_full_poly()`	length → model (9 / 19); raises otherwise
`soloff_velocity()`	3C reconstruction on the PIV grid
`make_world_grid()` / `_stereo_recon_grid()`	overlap world grid (full / PIV resolution)
`backproject_image()` / `backproject_folder()`	dewarp raw images to world coordinates
`compute_disparity()` / `disparity_to_delta_z()`	disparity-based \(\Delta z\) self-calibration (§10)

References¶

S. M. Soloff, R. J. Adrian, Z.-C. Liu (1997), “Distortion compensation for generalized stereoscopic particle image velocimetry,” Meas. Sci. Technol. 8, 1441–1454. — The mapping function \(M(X, Y, Z)\) and its Jacobian \(\partial M/\partial Z\).
C. E. Willert (1997), “Stereoscopic digital particle image velocimetry for application in wind tunnel flows,” Meas. Sci. Technol. 8, 1465–1479. — Measuring disparity by dewarp-and-correlate.
B. Wieneke (2005), “Stereo-PIV using self-calibration on particle images,” Exp. Fluids 39, 267–280. — Interpreting the disparity as a \(\Delta z\) offset between the laser sheet and the calibration plane.