# 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](#stereo-sec-4)). The model is selected by
the `mode` argument of {func}`~dpivsoft.stereo_DPIV.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:
{func}`~dpivsoft.stereo_DPIV._poly9_design`; fit:
{func}`~dpivsoft.stereo_DPIV._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{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}
$$

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:
{func}`~dpivsoft.stereo_DPIV._poly19_design`; fit:
{func}`~dpivsoft.stereo_DPIV._fit_poly19`. The GUI refuses to run the full
model with fewer than 3 distinct $Z$ planes per camera.

---

(stereo-sec-4)=
## 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
  ({func}`~dpivsoft.stereo_DPIV.make_world_grid`).
- **`map_*`** — world → pixel (the **Soloff map**):
  $(x_\mathrm{px}, y_\mathrm{px}) = \mathrm{poly}(X, Y, Z)$.
  Used by {func}`~dpivsoft.stereo_DPIV.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.

---

(stereo-sec-5)=
## 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{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}
$$

**Full (19-term), at any $Z$:**

$$
\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}
$$

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](#stereo-sec-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:

1. Map the grid point into each camera's raw image
   ({func}`~dpivsoft.stereo_DPIV._poly_eval` on `map_*`) and interpolate that
   camera's measured PIV displacement there (`griddata`).
2. Build each camera's Jacobian
   ({func}`~dpivsoft.stereo_DPIV._poly_jacobian` on `map_*`).

For camera $i$, the Jacobian links world displacement to that camera's measured
image displacement:

$$
\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}
$$

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

$$
\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}
$$

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
{func}`~dpivsoft.stereo_DPIV._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.

---

(stereo-sec-9)=
## 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`.
- **{func}`~dpivsoft.stereo_DPIV.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.

---

(stereo-sec-10)=
## 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](#stereo-sec-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**.

{func}`~dpivsoft.stereo_DPIV.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](#stereo-sec-5)) into the
in-plane $2 \times 2$ block and the parallax column:

$$
\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)
$$

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,
{func}`~dpivsoft.stereo_DPIV.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 |
|----------|------|
| {func}`~dpivsoft.stereo_DPIV.stereo_calibration` | fit both cameras, both directions; returns the 8 coefficient vectors and a results string |
| {func}`~dpivsoft.stereo_DPIV._fit_poly9` / {func}`~dpivsoft.stereo_DPIV._fit_poly19` | least-squares fit of one map (reduced / full) |
| {func}`~dpivsoft.stereo_DPIV._poly9_design` / {func}`~dpivsoft.stereo_DPIV._poly19_design` | design matrices (also used for the fit RMS) |
| {func}`~dpivsoft.stereo_DPIV._poly9_eval` / {func}`~dpivsoft.stereo_DPIV._poly19_eval` / {func}`~dpivsoft.stereo_DPIV._poly_eval` | evaluate a map (reduced / full / dispatch) |
| {func}`~dpivsoft.stereo_DPIV._poly9_jacobian` / {func}`~dpivsoft.stereo_DPIV._poly19_jacobian` / {func}`~dpivsoft.stereo_DPIV._poly_jacobian` | map derivatives (reduced / full / dispatch) |
| {func}`~dpivsoft.stereo_DPIV._is_full_poly` | length → model (9 / 19); raises otherwise |
| {func}`~dpivsoft.stereo_DPIV.soloff_velocity` | 3C reconstruction on the PIV grid |
| {func}`~dpivsoft.stereo_DPIV.make_world_grid` / {func}`~dpivsoft.stereo_DPIV._stereo_recon_grid` | overlap world grid (full / PIV resolution) |
| {func}`~dpivsoft.stereo_DPIV.backproject_image` / {func}`~dpivsoft.stereo_DPIV.backproject_folder` | dewarp raw images to world coordinates |
| {func}`~dpivsoft.stereo_DPIV.compute_disparity` / {func}`~dpivsoft.stereo_DPIV.disparity_to_delta_z` | disparity-based $\Delta z$ self-calibration ([§10](#stereo-sec-10)) |

---

## References

1. S. M. Soloff, R. J. Adrian, Z.-C. Liu (1997), "[Distortion compensation for
   generalized stereoscopic particle image
   velocimetry](https://doi.org/10.1088/0957-0233/8/12/008)," *Meas. Sci.
   Technol.* **8**, 1441–1454. — The mapping function $M(X, Y, Z)$ and its
   Jacobian $\partial M/\partial Z$.
2. C. E. Willert (1997), "[Stereoscopic digital particle image velocimetry for
   application in wind tunnel
   flows](https://doi.org/10.1088/0957-0233/8/12/010)," *Meas. Sci. Technol.*
   **8**, 1465–1479. — Measuring disparity by dewarp-and-correlate.
3. B. Wieneke (2005), "[Stereo-PIV using self-calibration on particle
   images](https://doi.org/10.1007/s00348-005-0962-z)," *Exp. Fluids* **39**,
   267–280. — Interpreting the disparity as a $\Delta z$ offset between the
   laser sheet and the calibration plane.