FloodTruth — LSPIV — stream discharge from video

1 · Method overview

Large-Scale Particle Image Velocimetry (LSPIV) adapts laboratory PIV — originally a two-flash, seeded-tracer technique — to natural rivers, where naturally-occurring surface tracers (foam, ripples, sun-glint patches, leaves, suspended debris) replace seeded particles. The technique was introduced by Fujita, Muste & Kruger (1998) and matured into routine practice through Muste et al. (2008), Le Coz et al. (2010), and free toolkits such as RIVeR (Patalano & García 2017), Fudaa-LSPIV (Le Coz et al. 2014) and KLT-IV (Perks 2020).

The pipeline implemented here is the standard four-stage LSPIV workflow: (i) camera-shake stabilisation against stationary bank features, (ii) dense optical flow on the water surface, (iii) time-averaging to suppress wave noise, (iv) velocity-area integration along a surveyed cross-section.

2 · Frame extraction & sampling

The video is decoded by stepping HTMLVideoElement.currentTime at fixed intervals $\Delta t = 1/f_s$ and copying each decoded frame into an offscreen canvas. The default sampling rate is $f_s = 10$ Hz, which is well above the ~1–3 Hz wave frequency of typical surface ripples and well below the ~30 Hz native rate of phone cameras.

3 · Image stabilisation (bank registration)

Hand-held footage shakes. To recover a clean velocity field on the water, every frame is geometrically registered to the first using stationary bank features:

1. Shi–Tomasi corner detector (Shi & Tomasi 1994) on the bank-mask region gives ~200 strong corner features:

$\displaystyle R(x,y) = \min\!\big(\lambda_1(x,y),\,\lambda_2(x,y)\big)$

where $\lambda_1, \lambda_2$ are the eigenvalues of the local image-gradient structure tensor.

2. Lucas–Kanade pyramidal optical flow (Lucas & Kanade 1981; Bouguet 2001) tracks each corner from frame $k$ to frame $k\!+\!1$, solving the brightness-constancy equation $I_x u + I_y v + I_t = 0$ in least-squares form over a local window:

$\displaystyle \begin{bmatrix} u\\v \end{bmatrix} = \big(A^\top A\big)^{-1} A^\top b,\qquad A=\begin{bmatrix} I_x & I_y \end{bmatrix},\quad b = -I_t$

3. Robust homography $H_{k}$ from the matched bank features via RANSAC, mapping frame $k$ back onto the reference frame (frame 0):

$\displaystyle \begin{bmatrix} x'\\y'\\1 \end{bmatrix} \;\sim\; H \begin{bmatrix} x\\y\\1 \end{bmatrix}$

Each frame is warped by $H_k$ before optical flow on the water is computed, removing camera translation and small rotations. (Strong perspective changes would require true orthorectification — see §8 limitations.)

4 · Dense optical flow on the water (Farnebäck)

Between each consecutive stabilised frame pair $(I_k, I_{k+1})$, the Farnebäck (2003) dense optical flow algorithm fits a local 2-D quadratic intensity model $I(\mathbf{x}) \approx \mathbf{x}^\top A \mathbf{x} + \mathbf{b}^\top \mathbf{x} + c$ inside a window centred on every pixel. Under a small displacement $\mathbf{d}$ the coefficients of the model in the next frame relate to those in the current frame via:

$\displaystyle A_2 \mathbf{d} = -\tfrac{1}{2}(\mathbf{b}_2 - \mathbf{b}_1)$

which is solved at every pixel, yielding a dense flow field $(u_{k}(x,y), v_{k}(x,y))$ in pixels per frame. The OpenCV implementation uses a 3-level pyramid, window size 15, 3 iterations, polynomial expansion size 5, and standard deviation 1.2.

Contrast gain. OpenCV's solve adds a small fixed regularising constant, so on a faint water texture the recovered displacement is biased toward zero — measured on a synthetic translating texture: −12% at a grey-level SD of 8, −0.3% at SD 20, 0% at SD 50. Before the flow step every stabilised frame gets the same affine intensity gain (target SD 40, measured on the first frame over the water, capped at ×8), which removes the scale dependence without violating brightness constancy. The gain applied is shown in the results. Frame pairs whose two frames are identical (a seek that returned the previous frame) are excluded from the average rather than counted as zero motion.

5 · Time-averaging

Wave-induced surface motion is approximately zero-mean over time, so averaging across all $N\!-\!1$ frame pairs suppresses wave noise and isolates the mean advective flow:

$\displaystyle \bar u(x,y) = \frac{1}{N-1}\sum_{k=1}^{N-1} u_k(x,y),\quad \bar v(x,y) = \frac{1}{N-1}\sum_{k=1}^{N-1} v_k(x,y)$

The pixel-frame flow is converted to physical surface velocity using the calibration:

$\displaystyle V_s(x,y) = \sqrt{\bar u^2 + \bar v^2}\;\cdot\; m_{\rm px}\;\cdot\; f_s \quad\text{[m/s]}$

where $m_{\rm px}$ is metres per pixel from §7 calibration and $f_s$ is the sampling rate.

6 · Pixel-to-metre calibration (nadir mode)

A known-length reference visible in the frame (a 1-m PVC pipe, surveyed bank-pin spacing) gives the image scale directly. Tapping its endpoints yields a pixel distance $L_{\rm px}$, so:

$\displaystyle m_{\rm px} = \dfrac{L_{\rm m}}{L_{\rm px}}\quad\text{[m/px]}$

This is correct only if the reference lies in the same plane as the water surface and the camera is roughly nadir-pointing. A 30°–60° tilt (typical bridge-railing shot) produces both a global cosine-error of $1-\cos\theta$ and a perspective stretch that varies across the frame — for those situations use ortho mode (§6b).

6b · Projective orthorectification (ortho mode)

When the camera is oblique to the water plane, every pixel sees a different ground footprint and a single $m_{\rm px}$ is meaningless. The standard photogrammetric fix is a plane-to-plane homography: a 3×3 matrix $\mathbf{H}$ relating homogeneous image coordinates $(u,v,1)^T$ to homogeneous world-plane coordinates $(X,Y,1)^T$:

$\displaystyle s\,\begin{bmatrix} X\\Y\\1 \end{bmatrix} = \mathbf{H}\,\begin{bmatrix} u\\v\\1 \end{bmatrix}\!,\qquad \mathbf{H}=\begin{bmatrix} h_{11} & h_{12} & h_{13}\\ h_{21} & h_{22} & h_{23}\\ h_{31} & h_{32} & 1 \end{bmatrix}$

Eight unknowns ($h_{33}$ is fixed to 1 to remove projective scale ambiguity), so four point correspondences are sufficient. The user taps four image points $(u_i,v_i)$ that they know form a rectangle of known width $W$ and length $L$ on the ground plane; the corresponding world-plane targets are the rectangle corners $(0,0)$, $(W,0)$, $(W,L)$, $(0,L)$.

The direct linear transform (DLT) writes the projection equation as two linear equations per correspondence in the 8 unknowns $\mathbf{h}=(h_{11},\dots,h_{32})^T$:

$\displaystyle \begin{pmatrix} u_i & v_i & 1 & 0 & 0 & 0 & -X_i u_i & -X_i v_i \\ 0 & 0 & 0 & u_i & v_i & 1 & -Y_i u_i & -Y_i v_i \end{pmatrix}\mathbf{h} = \begin{pmatrix} X_i\\ Y_i \end{pmatrix}$

Stacking four correspondences gives an 8×8 linear system $A\mathbf{h}=\mathbf{b}$ which is solved here by Gauss–Jordan elimination with partial pivoting. The system is well-conditioned for any non-degenerate (i.e. non-collinear, non-coincident) point quadruple. Once $\mathbf{H}$ is known, every video frame is warped to the rectified ground-plane view via OpenCV.js cv.warpPerspective, and the existing Farnebäck pipeline runs on the rectified frames — every pixel is now metric.

For N ≥ 4 surveyed GCPs (table mode, v0.4) the same DLT is solved as an over-determined system $A\mathbf{h}=\mathbf{b}$ of size $2N\times 8$, in the least-squares sense via the normal equations $A^{\!\top}\!A\,\mathbf{h}=A^{\!\top}\mathbf{b}$. Hartley's isotropic normalisation is applied to both the image and world point sets first — translate each to its centroid and scale so the RMS distance to the centroid is $\sqrt{2}$ — which keeps the conditioning of $A^{\!\top}\!A$ tight even when the world coordinates are large (e.g., UTM eastings in the $10^5$–$10^6$ m range). The recovered $\mathbf{H}$ is denormalised by $\mathbf{H} = T_{\!w}^{-1}\,\hat{\mathbf{H}}\,T_{\!i}$, and the per-GCP reprojection residuals (image-pixel distance between the tapped $(u_i,v_i)$ and the back-projected surveyed $(X_i,Y_i)$) are reported in the table for blunder detection. Reference: Hartley & Zisserman (2004) §4.4.4 "Normalising transformations."

The rectified output has user-chosen width $W_{\rm out}$ pixels and height $H_{\rm out} = W_{\rm out}\cdot L/W$ pixels (preserving the rectangle's aspect ratio), so the pixel-to-metre scale on the rectified frame is exactly:

$\displaystyle m_{\rm px}^{\rm ortho} = \dfrac{W}{W_{\rm out}} = \dfrac{L}{H_{\rm out}}\quad\text{[m/px]}$

Manual two-point calibration is unnecessary in ortho mode, and the cross-section line is tapped on the rectified preview, so its endpoints are already in metric ground coordinates. Reference: Hartley & Zisserman (2004) §4 "Estimation — 2D projective transformations," in Multiple View Geometry in Computer Vision (2nd ed.), Cambridge.

• Planarity assumption. The homography rectifies a single plane only — the water surface plane. Any feature off that plane (the rectangle had better be on the water plane, or directly above/below it on a flat bank) projects onto the wrong ground location and induces parallax error. Seasonal stage change between when the GCP rectangle was surveyed and when the video was shot does the same.
• Wave non-planarity. Orthorectification only corrects in-plane perspective. It cannot recover the 3-D shape of a non-planar instantaneous water surface — wave crests are still seen obliquely after the warp, they just appear in the wrong ground-plane location for an instant. This is acceptable in LSPIV because the time-average over the full clip smooths surface waves out, and we only report the time-mean velocity field.
• Resampling blur. Bilinear warping smooths the frame slightly. With sensible output sizes ($W_{\rm out}$ ≈ short-side native resolution) this is negligible compared to optical-flow noise.

7 · Cross-section velocity-area integration

Discharge is computed by the USGS mid-section method (Rantz 1982) along the surveyed cross-section. The user-tapped endpoints define the line; depths $d_i$ at stations $s_i$ come from bathymetric survey. The time-averaged surface speed $V_{s,i}$ at each station is sampled by bilinear interpolation of the $V_s(x,y)$ field at the pixel coordinates corresponding to $(s_i)$ on the cross-section line.

Surface velocity is converted to depth-averaged velocity using the velocity coefficient $\alpha$:

$\displaystyle \bar V_i = \alpha\,V_{s,i}$

The mid-section subwidth $\Delta w_i$ for each interior station is half the centre-to-centre distance to its neighbours; for the two end stations (which are at the waterline with $d=0$, $V=0$) the width is the half-distance to the single adjacent interior station. Per-segment discharge:

$\displaystyle q_i = \alpha\,V_{s,i}\,d_i\,\Delta w_i$

and total discharge:

$\displaystyle Q = \sum_i q_i$

8 · Known biases & limitations

• Cosine error from oblique angle (nadir mode only). An oblique camera projects pixels at the water surface onto a tilted plane, so $m_{\rm px}$ varies across the frame and the flat-field calibration is wrong. A 20° tilt produces ~6% bias in nadir regions and much larger errors near the far bank. Switch to ortho mode (§6b, v0.2) when the camera is at typical bridge-railing angles (30°–60° below horizontal).
• Surface-to-depth-averaged coefficient ($\alpha$). The standard value 0.85 (Rantz 1982; USGS Water-Supply Paper 2175) is a regression mean. Real $\alpha$ varies with bed roughness, depth, and channel shape: laboratory and field studies show 0.7–0.95. Hauet et al. (2018) found systematic 5–15% discharge bias at $\alpha=0.85$ relative to ADCP reference at 67 USGS sites.
• Dead-water near banks. Optical flow reports zero or noisy velocities along the channel margins where the water is shallow and there are no tracers. The mid-section integration handles this correctly only if the user enters station-0 and station-end as $d=0$ on the depth profile.
• Low tracer density. Glassy, smooth-water reaches with no foam, ripples, or debris produce flat optical-flow output (the gradient constraint is degenerate). Throw a handful of biodegradable tracers (sawdust, leaves) upstream before recording, or accept that LSPIV is not the right method for that reach.
• Shake registration assumes rigid bank. If the entire frame is moving water, the homography fit will lock onto the dominant flow and erase the very signal you want. Mitigation: increase bank-mask exclusion, ensure visible stationary features (rocks, vegetation, banks) in the top + bottom strips.
• Wave-only motion. On wind-rippled lakes or backwater, surface tracers move with the wave orbital velocity, not with bulk transport. Net advection averages to zero over enough frames, but the SNR is poor. LSPIV is most reliable in moving rivers.
• Sun glint / reflections. Specular reflections from sun, sky, and overhanging vegetation create bright features that are stationary in earth coordinates but move when the camera tilts — the homography may misregister against them. Polarising filter or overcast conditions help.
• Aliasing at high velocity (auto-detected, v0.2). Farnebäck assumes small per-frame displacements — <15 px is the routine working range, beyond ~25 px the polynomial-expansion residuals saturate and the reported flow is biased low and noisy (see §3 of the existing methodology, citing Farnebäck 2003). After time-averaging the flow field, the app computes the 95th-percentile of $|\mathbf{v}|$ in pixels-per-frame across the analysed region (banks excluded). If >15 px/frame the user is warned that velocities are biased low and a recommended sampling rate $f_{\rm rec} = f_s\cdot p_{95}/12$ is reported. If >25 px/frame the result is flagged "unreliable" (red) and $Q$ is tagged biased-low. Mitigation: increase target fps (try 20–30 Hz), get closer to the water so $m_{\rm px}$ is larger per pixel, or both.
• Compressed video artefacts. Heavy H.264 compression on phone video erases fine-scale tracers and creates block-edge "ghost flow." Record at the highest bitrate / least compression setting available.

9 · References

• Fujita, I., Muste, M., Kruger, A. (1998). Large-scale particle image velocimetry for flow analysis in hydraulic engineering applications. Journal of Hydraulic Research 36(3), 397–414.
• Muste, M., Fujita, I., Hauet, A. (2008). Large-scale particle image velocimetry for measurements in riverine environments. Water Resources Research 44, W00D19.
• Le Coz, J., Hauet, A., Pierrefeu, G., Dramais, G., Camenen, B. (2010). Performance of image-based velocimetry (LSPIV) applied to flash-flood discharge measurements in Mediterranean rivers. Journal of Hydrology 394, 42–52.
• Patalano, A., García, C.M., Rodríguez, A. (2017). Rectification of image velocity results (RIVeR): a simple and user-friendly toolbox for large-scale water surface PIV and PTV. Computers & Geosciences 109, 323–330.
• Perks, M.T. (2020). KLT-IV v1.0: image velocimetry software for use with fixed and mobile platforms. Geoscientific Model Development 13, 6111–6130.
• Hauet, A., Morlot, T., Daubagnan, L. (2018). Velocity profile and depth-averaged to surface velocity in natural streams: a review over a large sample of rivers. E3S Web of Conferences 40, 06015.
• Farnebäck, G. (2003). Two-frame motion estimation based on polynomial expansion. SCIA 2003: Image Analysis, LNCS 2749, 363–370.
• Lucas, B.D., Kanade, T. (1981). An iterative image registration technique with an application to stereo vision. Proc. IJCAI, 674–679.
• Shi, J., Tomasi, C. (1994). Good features to track. Proc. IEEE CVPR, 593–600.
• Bouguet, J.-Y. (2001). Pyramidal implementation of the affine Lucas-Kanade feature tracker. Intel Corp. Technical Report.
• Rantz, S.E. (1982). Measurement and computation of streamflow. Vol. 1: measurement of stage and discharge. USGS Water-Supply Paper 2175.
• Hartley, R., Zisserman, A. (2004). Multiple View Geometry in Computer Vision (2nd ed.), Cambridge University Press, §4 — 2D projective transformations and the DLT algorithm.
• Abdel-Aziz, Y.I., Karara, H.M. (1971). Direct linear transformation from comparator coordinates into object space coordinates in close-range photogrammetry. Proc. ASP/UI Symposium on Close-Range Photogrammetry, 1–18.