O R I G I N A L S
D. Douxchamps Æ D. Devriendt Æ H. Capart Æ C. Craeye B. Macq Æ Y. Zech
Stereoscopic and velocimetric reconstructions of the free surface
topography of antidune flows
Received: 14 February 2004 / Revised: 10 March 2005 / Accepted: 17 March 2005 / Published online: 15 June 2005 Springer-Verlag 2005
Abstract Imaging methods developed to characterise the oscillatory free surface of rapid flows are presented and applied to torrential currents over sediment antidunes. The aim is to obtain high-resolution relief maps of the free surface topography. Two measurement principles are tested, both based on the imaging of floating tracers dispersed on the rapidly flowing surface. The first tech-nique involves direct stereoscopic measurements. The second technique is indirect, and exploits a Bernoulli relation to derive surface elevations from the horizontal velocity field acquired using a single camera. Special attention is paid to error estimation and control. Relief maps obtained for various bedform patterns are pre-sented, allowing comparison between the two tech-niques.
1 Introduction
Free surface measurements constitute an important challenge in experimental fluid mechanics. The free surface is often of interest because of the key role it plays in the flow phenomena under consideration (e.g. Ham-mack et al. 1989; Ja¨hne et al. 1994; Lang and Gharib
2000). Moreover, it may be the only flow feature, which
is readily accessible to measurements by imaging or re-mote sensing methods (e.g. Stilwell and Pilon 1974; Nicolas et al. 1997; Craeye et al. 1999). Both circum-stances are encountered in the present study of antidune flows.
Antidunes are trains of bed waves (Gilbert 1914; Kennedy 1963; Reynolds 1965; Allen 1968) which ap-pear when rapid, shallow currents flow over coarse granular material, and are characterised by in-phase coupling between the oscillatory sediment bed and the water free surface. As shown in Fig.1, they can be ob-served both in the field (Alexander and Fielding 1997; Blair 2000) and in the laboratory (Middleton 1965). Along with other types of fluvial bedforms, they are of interest to geomorphologists and hydraulic engineers (Shaw and Kellerhals 1977). Antidunes are notable in particular for their evanescent character: they vanish rapidly once the flow wanes, and leave few lasting traces aside from bedding and grain sorting effects. As a result, their geometrical configuration is best studied when the flow is active.
Long-crested or short-crested, arranged in regular arrays or in narrow trains of peaks and troughs, antid-unes come in a variety of patterns. In the present study, a characterisation of such patterns is sought through measurements of the water free surface topography. Reflection and refraction by the rough free surface im-pedes visual access to the underlying sediment bed, hence no direct measurement of the bottom topography is possible. Since the two surfaces are locked in phase with each other, however, the water surface topography provides an indirect image of the bedform pattern.
Measurements of water surface topography most often involve point sensors, placed in multi-sensor arrays at fixed locations or scanned across the surface (e.g. Hammack et al.1989; Wessels et al.1989). Sensors used to measure water elevation include resistive or capacitive gauges, pressure transducers and acoustic beams. Resistive gauges are inapplicable in the present case because of the high sensitivity of antidune flows to intruding objects. More generally, the spatial and
D. Douxchamps Æ C. Craeye Æ B. Macq
Communication and Remote Sensing Laboratory, Universite´ catholique de Louvain, 3, place du Levant, 1348 Louvain-la-Neuve, Belgium
D. Devriendt Æ Y. Zech
Department of Civil and Environmental Engineering, Universite´ catholique de Louvain, 1, place du Levant, 1348 Louvain-la-Neuve, Belgium
H. Capart (&)
Department of Civil Engineering and Hydrotech Research Institute, National Taiwan University, 10617 Taipei, Taiwan
E-mail: hcapart@yahoo.com DOI 10.1007/s00348-005-0983-7
temporal resolution of point sensors is limited by the number of available devices or the time required to scan a sensor across the field of interest. This makes them unsuitable for the transient (on a time scale of a few seconds), spatially varied surfaces of flows over evolving bedforms.
Because they permit non-intrusive, fast whole-field measurements, imaging techniques constitute attractive alternatives to point gauges and sensors. Described in a review by Ja¨hne et al. (1994), a variety of imaging principles have been proposed to characterise water surface shapes. Based on photogrammetric techniques widely applied to land topography, one approach is to measure water surface heights using stereo photography (Banner et al. 1989; Shemdin and Tran 1992). This re-quires corresponding features or regions on the surface to be matched between distinct views. For water sur-faces, this correspondence problem is difficult to solve reliably because of the specular nature of light reflection at the free surface (Ja¨hne et al.1994).
For this reason, most investigators have turned to measuring surface slope rather than surface height. Approaches based on shading, reflection and refraction have been developed for this purpose, and are reviewed in Ja¨hne et al. (1994). Recent measurements based on reflection are documented in Craeye et al. (1999) and Dabiri and Gharib (2001), while recent studies relying on refraction (shadowgraph) include Weigand (1996) and Lang and Gharib (2000). These techniques,
however, are subject to various limitations. Shape-from-shading does not work for transparent, specularly reflecting surfaces. On the other hand, shape-from-reflection is restricted to a narrow slope range. Finally, shape-from-refraction requires light rays to pass through the top and the bottom of the water layer. These limi-tations make the techniques inapplicable to the present water flows featuring rough and highly varied free sur-faces, and propagating over opaque sediment beds.
Like the earlier stereo techniques, the approach ex-plored in the present work relies on the matching and tracking of surface features. To address the correspon-dence problem, floating particles are dispersed on the water free surface, furnishing point-like features that look the same on distinct frames. Images of these floating tracers are then exploited using two different reconstruction principles. The first is a classical stereo-scopic principle based on the matching of particles on image pairs acquired from two cameras. The second is an original velocimetric principle requiring only a single vertical camera: the surface velocity field is first acquired by particle tracking velocimetry, then converted into a vertical elevation map using a Bernoulli relation derived from the fluid mechanics of the water free surface.
The parallel development of two distinct techniques was motivated by the following considerations. First, it makes it possible to weigh the respective advantages of the methods and evaluate their applicability to more complex situations. Second, an assessment of the valid-ity of both techniques can be obtained by comparing their results. To permit cross-validation, the two meth-ods are applied to the same experiments, while each re-lies on its own separate camera footage and data analysis pipeline. Preliminary results of the present re-search effort were reported in Devriendt et al. (1998) and Douxchamps et al. (2000).
The paper is organised as follows. In Sect. 2, the experimental setup and camera systems used for the measurements are presented. Then, basic particle imag-ing algorithms used as buildimag-ing blocks by the two reconstruction methods are outlined. Sections 3 and 4
are devoted to a detailed presentation of each of the two reconstruction methods. This includes special proce-dures for the estimation and filtering of measurement errors at each step of the analysis. For the stereo prin-ciple, the surface reconstruction and error estimation procedures are validated using a solid surface of known shape. Finally, the methods are applied to the free sur-faces of water flows over antidunes. Results from the stereo and velocity-based methods are shown, compared and discussed in the last section preceding the overall conclusions of the work.
2 Experimental setup and imaging configuration
Experiments are conducted in a hydraulic flume having the following dimensions: length=6 m; width=50 cm; side-wall height=40 cm. The flume is tilted to obtain a
Fig. 1 Antidune flows on a beach of Eastern Taiwan (top) and in the Louvain laboratory flume (bottom). Photographs by B. Spinewine and H. Capart
bottom slope of 1%. A 5-cm deep layer of loose sedi-ment covers the flume bottom and is replenished during the flow by an upstream silo. Sand of nearly uniform size distribution is used as sediment material, and has the following characteristics: mean grain diame-ter=1.65 mm; density=2,615 kg/m3. For the proof-of-concept experiments presented hereafter, the water inflow was not tightly controlled but rather loosely operated in order to produce a variety of patterns. Water discharges in the flume ranged from 15 l/s to 35 l/ s, yielding flow depths of 3–6 cm and Froude numbers in the range Fr=1.4–2. For such high Froude numbers, the water free surface responds strongly to the underlying oscillatory topography of the bedforms, up to the point of breaking at the wave crests when antidunes are fully developed. The experimental parameters of the three runs selected for analysis are listed in Table 1.
Starting from a plane bed, the antidunes emerge as longitudinal trains of crests and troughs initiating from downstream but stationary in phase. The antidunes are observed to respond to transient changes in the flow rate (both increase and decrease) by temporarily growing in amplitude. Amplitude responses and gradual shifts in pattern occur on a time scale of tens of seconds. By contrast, on the shorter time scale corresponding to the image acquisition (of the order of 2 s), the hydrody-namics can be assumed to be quasi-steady and this is exploited hereafter to derive a single surface from each measurement sequence. The flow is however observed to exhibit small unsteady pulses, and this physical source of noise may slightly affect the results concurrently with measurement errors.
The short wavelength features of the water surface itself are not suitable for either matching or tracking. Specular reflection changes the appearance of these features when viewed from different angles. They also continuously evolve in time under the action of capil-larity and gravity. Floating particles dispersed on the water free surface are used instead as tracers. The tracer particles are white wooden pearls 9 mm in diameter. The measurement section is placed some 2 m upstream of the flume outlet. Further upstream and moments before image acquisition, a uniform dispersion of floats is dropped onto the mean surface by means of staggered metal meshes.
Image sequences are obtained using digital cameras placed above the flow. The velocimetric and stereoscopic methods require rather different image acquisition sys-tems, hence some care is necessary in order to operate them simultaneously and observe the same scene. Two commercial digital cameras (miniDV, PAL, 25 frames per second) are used for the stereoscopic measurements. These cameras offer good image resolution (768/ 576 pixels) but cannot be synchronised with each other during the acquisition. To avoid motion-stereo ambi-guity, it will thus be necessary to synchronise them a posteriori using an interpolation procedure (see below). The stereo cameras are placed above the flow with ob-lique optical axes contained in a vertical plane parallel to the direction of flow (see Fig. 2). For the velocimetric measurements, on the other hand, a high frame rate camera (250 frames per second; resolution of 239/ 192 pixels) is placed directly over the flow with a nearly vertical optical axis. Due to the high frame rate, this
Table 1 Experimental parameters and error estimates
Run no. 0 3 5 6
Surface pattern Sinusoidal plate Rolls Zig-zag Narrow train
Mean surface velocity (mm/s) NA 970 1,180 980
Mean flow depth (mm) NA 40 52 31
Wavelength (mm) 41.1 320 410 330
Surface elevation range gmax gmin(mm) 7.0 40 40 40
Stereoscopic method
Number of matched particles pairs 4,725 13,300 11,100 11,500
Point elevation rms error ^rg (mm) 0.9 2.8 2.8 2.8
Ratio of correct matches P(x1) (%) 93 58 57 61
Noise error attenuation error (mm) 0.27 0.51 0.52 0.50
Reconstruction rms error (mm) 0.75 0.96 1.07 0.99
Overall elevation rms error (mm) 0.84 1.2 1.3 1.2
Velocimetric method
Number of tracked particle pairs NA 48,600 38,700 47,300
Point position rms error rx(mm) NA 0.31 0.31 0.29
Noise error attenuation error (mm/s) (trajectory filtering) NA 8.4 8.4 7.9 Noise error attenuation error (mm/s) (transverse averaging) NA 7.9 8.7 8.4
Overall velocity rms error (mm/s) NA 16.3 17.1 16.3
Velocity error converted to elevation (mm) NA 1.6 2.1 1.6
Second order rms perturbation (mm) NA 2.1 1.7 1.7
Overall elevation rms error (mm) NA 2.7 2.7 2.4
Comparison
Predicted rms discrepancy (mm) 0.84 2.9 3.0 2.7
Measured rms discrepancy (mm) 1.11 3.8 5.1 4.7
camera requires strong lighting, obtained with four 2 kW light sources. Such powerful lighting saturates the commercial cameras even at maximum shutter speed, and these have to be fitted with dimming filters.
After positioning, the viewpoints of all three cameras are determined by placing a calibration target in their shared viewing volume. This is essential for stereo reconstruction and allows the results of the two methods to be obtained in the same three-dimensional referential. For ease of reference in the next sections, we adopt la-bels V=L, R and T to identify the three viewpoints. Labels L and R designate the left and right oblique cameras used for stereoscopic reconstruction, while label T refers to the top camera used for velocimetric recon-struction. The camera configuration shown in Fig.2b is reconstructed based on actual calibrated viewpoints.
3 Particle imaging algorithms
While the two reconstruction methods are distinct, they rely on shared building blocks provided by general-purpose particle imaging algorithms. These algorithms achieve three basic tasks: (1) identification of particles on the two-dimensional images; (2) tracking of particles in the image plane; (3) establishment of a geometrical correspondence between 2D image coordinates and positions within the 3D viewing volume. The basic algorithms we use for these purposes are either standard
or well-documented elsewhere, and are only summarised in the present section.
3.1 Particle positioning
Consider image sequences acquired under any of the three viewpoints V=L, R or T. For both reconstruction techniques, the first step of the analysis is the localisation of particle centres on individual images. One seeks to extract for each viewpoint V a set of positions fRðVÞi;mg where Rði;mV)¼ ðXi;mðV); Yi;mðV)Þ denotes the position of the i-th particle visible on i-the m-i-th image of i-the sequence. The white particles show up on the images as bright spots against a dark, relatively noisy background (see Fig.2). Such bright spots are first highlighted by convoluting the image with a laplacian-of-gaussian filter (Ja¨hne 1995). Local maxima of the highlighted images are then sought. The position of each maximum is finally refined to sub-pixel accuracy by way of a quadratic interpolation surface (for more details, see Capart et al. 2002). Illus-trating these steps, a filtered image with detected particle positions is shown in Fig.3a. The expected root mean square accuracy on the X and Y image coordinates obtained with such a procedure is of the order of 0.1–0.25 pixels (Veber et al.1997; Capart et al.2002).
3.2 Particle tracking
Particle tracking is required in the present work for both reconstruction methods. For velocimetric reconstruc-tion, particle tracking is used to acquire the horizontal
velocity field and constitutes an essential component of the method. For stereoscopic reconstruction, tracking is not an essential component, but is required in the present work because of the use of non-synchronised
cameras. A posteriori synchronisation is needed to avoid stereo-motion ambiguity, and is achieved by interpo-lating particle positions along their tracked trajectories. The particle-tracking step aims at establishing corre-spondence between particles seen on one image frame and particles seen on the next. Since particles are iden-tical, the information available is limited to the two sets of particle positions {RðVÞj;m} and {RðVÞj;mþ1} obtained for the two successive frames m and m+1. The tracking prob-lem is then one of finding the pairing ji) associating each particle j on frame m+1 to its most reasonable coun-terpart i on the previous frame m.
In the present work, this operation is performed using the Voronoı¨ algorithm of Capart et al. (2002). This algorithm tracks local patterns of neighbouring particles using Voronoı¨ cells as match templates. Voronoı¨ cells (see e.g. Okabe et al. 1992) are polygons surrounding each particle, encompassing the region of the plane closer to that particle than to any other. Qualitatively, the Voronoı¨ match algorithm pairs particles, which are characterised on two successive images by similar Vor-onoı¨ cells. The procedure is shown in Fig.3: Voronoı¨ diagrams constructed on two successive frames are shown on panel 3b and the corresponding displacement vectors are plotted on panel 3c. As shown in Fig.3d, displacements can finally be concatenated to extract particle trajectories over the entire image sequence. For certain image pairs, a small number of mismatches are obtained. They are filtered out using the procedures detailed in Capart et al. (2002).
3.3 Imaging geometry
One must now be able to relate the image coordinates RðV) ¼ ðXðV); YðV)Þ associated with each camera view-point V to a system of 3D world coordinates measured in a shared reference frame. The required transforma-tion can be modelled as a perspective projectransforma-tion (see Faugeras and Luong 2001). This amounts to specifying for each viewpoint V a matrix A(V) and a vector b(V)allowing image coordinates (X(V), Y(V)) to be com-puted from a known position according to
a XðVÞ YðVÞ 1 0 @ 1 A ¼ AðVÞ xy z 0 @ 1 A þ bðVÞ: ð1Þ
Conversely, a point having image coordinates (X(V), Y(V)) under viewpoint V belongs to a 3D ray defined by parametric equation
Fig. 3 Particle positioning and tracking operations applied to images from the top camera: a filtered image and particle positions (+); b Voronoı¨ tesselations constructed on particle positions identified on one frame (thin lines) and the next (thick lines); c displacement vectors obtained by matching the shapes of the Voronoı¨ cells (see Capart et al. 2002); d particle trajectories obtained by concatenating displacements over the whole sequence
rðVÞðaÞ ¼ pðVÞþ a qðVÞ; ð2Þ where scalar a is the free parameter. Vectors p(V)and q(V) denote the position of the projection centre and the ray direction, respectively, and are given by
pðVÞ¼ ðAðVÞÞ1bðVÞ qðVÞ¼ ðAðVÞÞ1 XðVÞ YðVÞ 1 0 @ 1 A: ð3a; bÞ A least squares camera calibration procedure is nee-ded to obtain matrix A(V) and vector b(V) for each viewpoint, as explained in more detail in Spinewine et al. (2003). In the present experiments, each of the view-points V=L, R and T is calibrated using some 200 points of known world and image coordinates. They are acquired from snapshots of a target placed within the viewing volume before or after the experimental runs.
4 Stereoscopic reconstruction
The first reconstruction method is a feature-based ste-reoscopic matching procedure used in digital photo-grammetry (see e.g. Wolf and Dewitt2000), and applied here to the floating tracer particles. The core of the method lies in the matching of rays associated with the particle images picked up by the two oblique cameras. This results in a three-dimensional cloud of points, from which the desired water free surface is recovered by surface fitting. In the following subsections, the basic ray intersection and matching operations are first outlined. Next, an original maximum-likelihood method devel-oped to filter out stereoscopic mismatches from the surface fit is presented. The procedures devised for error estimation and control are then described. To validate the methods, tests conducted with a surface of known shape are documented in the last subsection. Results for
the antidune flows will be presented further below after the velocimetric method has been described.
4.1 Ray intersection and matching
Consider two candidates for a stereoscopic match, characterised by coordinates Ri
(L)
and Rj (R)
on images seen by the left and right cameras (V=L, R). The cor-responding rays are given by parametric equations (see Eq. 2)
rðLÞi ðaÞ ¼ pðLÞi þ aqðLÞi rðRÞj ðbÞ ¼ pðRÞj þ bqðRÞj :
ð4a; bÞ With reference to Fig. 4, the points on the left and right rays corresponding to the minimum distance be-tween the two rays are parameterised by the values a, b, which are solutions to the system
qðiL) qðiL) qðiL) qðjR) qðjR) qðiL) qðjR) qðjR) 0 @ 1 A a b ¼ q ðLÞ i ðp ðRÞ j p ðLÞ i Þ qðRÞj ðpðRÞj pðLÞi Þ ! ; ð5Þ
where the dot denotes the usual dot product. This expresses the condition that the line segment linking the two points must be orthogonal to both rays. The mid-point ri,j and distance di,j between the two points of
closest encounter are then given by
ri;j¼ 1 2 r ðLÞ i ðaÞ þ r ðRÞ j ðbÞ di;j¼ rðRÞj ðbÞ r ðLÞ i ðaÞ : ð6a; bÞ
Due to the limited accuracy of the camera calibration and image plane measurements, rays corresponding to
Fig. 4 Intersection of two rays associated with the left and right camera viewpoints
the same physical point will not perfectly intersect. Midpoint ri,j constitutes an approximate intersection,
while the distance of closest encounter di,j provides an
indicator of the quality of the approximation. A procedure to select stereo pairs follows. Let {Ri
(L)
} and {Rj
(R)
} designate the two sets of particle positions picked up at the same time t by the left and right cam-eras. The task is to find the pairing ji) most likely to match together rays corresponding to one and the same physical particle. The stereo matching problem is in this regard identical to the velocimetric tracking task sket-ched earlier in Sect.3.2and documented in greater detail in Capart et al. (2002). It can be cast as an optimisation problem, whereby one seeks to minimise the objective function
X
i
di;jðiÞ ð7Þ
under the constraint that a given particle image can participate in only one stereo pair. This objective graph optimisation problem can be solved approximately using the Vogel algorithm: consider for each ray the best and second best match, then construct a reasonable global optimum by picking ray pairs in the order of maximum difference between first and second best choices. The set of 3D positions is finally obtained from (6a) as { ri,j(i)}.
The matching and intersection of rays associated with particle images thus yield the world coordinates of the physical particles.
Because the two cameras used for the antidune experiments could not be synchronised during acquisi-tion, the procedure above is not applied directly to the raw image positions. Instead, the particles seen under one of the two viewpoints are first tracked in the image plane using the methods of Sect. 3.2. Interpolation in time is then used to resample their positions at the in-stants of frame acquisition for the second viewpoint. This is needed to restore precise synchronicity of the position sets, and avoid the stereo-motion ambiguity, which would otherwise plague the results (Cameron
1952).
Pattern-based matching of stereo pairs similar to the Voronoı¨ cell tracking approach was not attempted. The reason lies in the distortion of cell shapes resulting from the widely different perspectives of the two stereo cam-eras. This problem could be solved by using a stereo setup with a shorter baseline or exploiting projective invariant characteristics of the point patterns. While this exceeds the scope of the present work, it could be a worthwhile avenue for further research.
4.2 Mismatch filtering using maximum-likelihood surface fitting
The above matching procedure results in sets of points ri,m positioned within the viewing volume at time
in-stants tm. Since in the present application the flow is
assumed to be quasi-steady over the duration of each
video acquisition, the different subsets ri,mcan further be
pooled into one single set rk. A continuous water free
surface elevation g(x, y) remains to be extracted from this cloud of points. If the matching process were per-fectly accurate, all points would belong to the desired surface and interpolation would suffice. In practice, however, various sources of error conspire to offset the points around the true surface, and some form of fitting is required.
In the present case, adequate fitting requires careful consideration of the structure of the measurement er-rors. The measured point set rk includes two distinct
subsets: (1) points resulting from correctly matched rays (i.e. they do indeed correspond to one and the same physical particle) but positioned with a limited accuracy; (2) points resulting from mismatched rays. It is reason-able to expect that points of the first subset will cluster around the physical surface with an error distribution that is approximately Gaussian, due both to position inaccuracy and to physical fluctuations. Mismatches, however, create a very different kind of background noise.
With reference to Fig.5a let us suppose for the sake of argument an idealised viewing geometry with focal
Fig. 5 Free surface reconstruction in the presence of stereo mismatches: a conceptual sketch; b comparison of mean (dashed line), median (thin line) and maximum-likelihood reconstructions (thick line) of the free surface based on stereo data points (dots)
points of the projection sent to infinity, and a uniform distribution of particles on a wavy surface. Geometri-cally, mismatches correspond to the establishment of a wrong correspondence between two rays that approxi-mately belong to the same epipolar plane (i.e. the plane containing the point of intersection of the two rays and the focal points of both projections). Among the family of such rays, each one is equally likely to produce a mismatch, hence mismatched points will be distributed randomly in the viewing volume rather than cluster around the physical surface. In an actual imaging geometry, the reasoning should still hold in an approx-imate way, hence the structure of the point distribution is expected to combine: (1) correct matches corrupted by Gaussian noise; (2) mismatches uniformly distributed in the viewing volume. In the presence of a high proportion of mismatches, highly biased results can be obtained if this structure is overlooked.
Rather than the whole set, let us consider the points belonging to a prismatic bin of dimensions Dx, Dy around position (x, y). Let x1 designate the correctly
matched state in which a particle can find itself, and x2
designate the mismatched state. P(x1) and P(x2) are the
corresponding probabilities. Based on the above model, the point elevations zk are randomly distributed with
probability density function given by
pðzjH1;H2Þ ¼ pðzjx1;H1ÞP ðx1Þ þ pðzjx2;H2ÞP ðx2Þ;
ð8Þ where H 1, H 2 designate the parameter vectors which
characterise the distribution of each class. This is a so-called mixture density, in which conditional densities p(z|x 1, H 1), p(z|x 2, H 2), are called component
densities and probabilities P(x 1), P(x2) are called the
mixing parameters (see e.g. Duda and Hart 1973). For the normal and uniform distributions postulated above, we have pðzjx1;H1Þ ¼ 1 rg ffiffiffiffiffiffi 2p p exp ðz gÞ 2 2r2 g ! pðzjx2;H2Þ ¼ ðzb zaÞ1; ð9a; bÞ
hence parameter vectors are given by H1= (g, rg) and
H 2= (za, zb), where g, rgare the mean and standard
deviations of the distribution of correct matches and za,
zb are the (possibly truncated) vertical limits of the
viewing volume at location (x, y). Symbol g has been used to designate the mean of the correct matches be-cause it designates the surface elevation that is ultimately sought.
The problem of finding the surface elevation has thus been cast into one of finding an optimal estimate for parameter g of a component density. Since it is not known a priori which points are correct matches and which are mismatches, the task at hand belongs to the category of unsupervised learning problems (Duda and Hart 1973). The elevation boundaries za and zb are
known and, by definition, P(x1)+P(x2)=1. Three
parameters thus remain to be determined: g, rg and
P(x1). In a manner analogous to the case of normal
mixtures treated by Duda and Hart (1973), maximum-likelihood estimates ^g;, ^rg and ^Pðx1Þ must satisfy
^ Pðx1Þ ¼ 1 n Xn k¼1 ^ Pðx1jzk; ^g; ^rgÞ; ð10aÞ ^ g¼ Pn k¼1 ^ Pðx1jzk; ^g; ^rgÞzk Pn k¼1 ^ Pðx1jzk; ^g; ^rgÞ ; ð10bÞ ^ r2g ¼ Pn k¼1 ^ Pðx1jzk; ^g; ^rgÞðzk ^gÞ2 Pn k¼1 ^ Pðx1jzk; ^g; ^rgÞ ; ð10cÞ
where, by Bayes’ rule (see Duda and Hart1973) ^
Pðx1jzk; ^g; ^rgÞ ¼
pðzkjx1; ^g; ^rgÞ^Pðx1Þ
pðzkjx1; ^g; ^rgÞ^Pðx1Þ þ pðzkjx2; za; zbÞ 1 ^Pðx1Þ
ð11Þ
and where pðzkjx1; ^g; ^rgÞ and p(zk|x2, za, zb) are given by
distributions (9). Despite the involved expressions, esti-mates (10a, b,c) can be interpreted rather simply as a frequency ratio, a sample mean and a sample variance, respectively, but with samples weighted according to the conditional probability that they belong to the class of correctly matched particles. A solution of the above set of implicit equations can be obtained by iterations. Initial estimates for ^g;^rg and ^Pðx1Þ can first be used to
evaluate ^Pðx1jzk; ^g; ^rgÞ using (11). Equations 10a, b,c
can then be used to update the estimates, and so on. While convergence to a single solution is not guar-anteed in all cases, the procedure was found to perform adequately in this application, provided reasonable ini-tial estimates are chosen. For the present antidune experiments, the converged estimates for ^rg and ^Pðx1Þ
are ^rg 3 mm and ^Pðx1Þ 0:6: The latter number is
especially significant: about 40% of the data points can be assessed to result from mismatches! This is due to the combination of limited image resolution with a dense dispersion of tracer particles on the free surface, creating many opportunities for epipolar ambiguity. The high likelihood of mismatches underscores the importance of filtering out their effect on the final measurements.
Figure5b presents results for a sample longitudinal slice of the measurements of run 3. Estimates obtained by binning the values, then taking the local mean or median rather than applying the above maximum-likelihood estimate are also shown. Although the median is slightly less affected, both these simpler estimates are seen to be significantly biased by the presence of uniformly dis-tributed mismatches in addition to the more classical Gaussian scatter. The maximum-likelihood Bayesian estimate exhibits a better behaviour, disregarding the
random mismatches and focusing on the densely clus-tered points. Similar techniques aimed at weighting down outlying measurements in the presence of non-Gaussian noise are used in high-energy physics (e.g. Fru¨hwirth et al.2000).
4.3 Balance of noise and attenuation errors
To carry out the above procedure, individual data points must be assigned to finite size bins. Overlapping square bins are chosen, with bin centres positioned on a regular grid xi·yj. These bins have a variable size, each one
encompassing a constant number of data points n. An optimal choice of this number n is based on the fol-lowing considerations. When applied to a spatially dis-tributed set of data corrupted by uncorrelated Gaussian noise, averaging produces two different effects: on the one hand, the Gaussian noise error diminishes in pro-portion to 1=pffiffiffin if only Gaussian error is present (i.e. there are no stereo mismatches); on the other hand, the coarser the averaging, the higher the attenuation of fine-grained details of the spatial field, inducing an attenua-tion error which grows with the bin size and number n of binned data. These two effects act at cross-purposes, and it is possible to choose a number n, which minimises the combined error.
Let ^gðnÞðx; yÞ designate the reconstructed surface resulting from assigning n data points to each bin and applying the Bayesian maximum-likelihood estimator of the previous section. An estimate of the residual Gaussian noise error is given by
ð^gðnÞ gðnÞÞ2 D E1=2 ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffir^g ^ Pðx1Þ n q ; ð12Þ
where the triangular brackets denote an ensemble aver-age, gðnÞis the averaged surface which would be obtained by bin-averaging in the absence of noise, ^rg is the
standard deviation of the Gaussian error as estimated from Eq. 10c, and ^Pðx1Þn is the number of Gaussian
distributed points estimated from Eq. 10a. Reduced estimate ^Pðx1Þn is used in Eq. 12 rather than the full
number of points n in order to disregard the non-Gaussian stereoscopic mismatches occurring with probability 1 P(x1).
An estimate of the attenuation error, on the other hand, can be obtained as follows. Let m designate the density of data points per unit surface obtained after collapsing all data into one single dataset. Then the likely bin size resulting from assigning a constant num-ber of points n to each bin is given by Dx¼ Dy ¼pffiffiffiffiffiffiffin=m. Assume now that the free surface is roughly harmonic, with variations in the x and y directions characterised by wavelengths kx ky [k] and amplitude a [g’], where
[k] and [g’] designate characteristic horizontal and ver-tical scales of the free surface oscillations. The averaging process can be idealised as a convolution of the har-monic surface with a square top-hat filter of size Dx ·
Dy. This results in attenuated amplitude aðnÞ given by (e.g. Ja¨hne1995) aðnÞ a ¼ sinc p Dx kx sinc pDy ky sinc2 p ffiffiffiffiffiffiffi n=m p ½k ! ; ð13Þ
where the sinc function is defined as sinc(x)=sin(x)/ x. Developing the sin c function by Taylor around the origin (up to second order), one can estimate the attenuation error due to the spatial averaging as
ðgðnÞ gÞ2 D E1=2 ¼ 12ða aðnÞÞ p 2 6 ½g0 m½k2n; ð14Þ
where the attenuation error is seen to grow in propor-tion to the number of binned points n. Finally, an approximate way of minimising the combined error is to set errors (12) and (14) equal to each other. For the conditions of the present antidune experiments (see Table1), this results in an optimal number of points n 50 and an average bin size Dx=Dy45 mm. The asso-ciated rms error on surface elevation due to either noise or attenuation is around 0.5 mm.
In addition to these two sources of error, which de-pend on the number of points n used for binning, the reconstruction generates a background error which is independent of n. For perfectly accurate elevation measurements (rg=0), the procedure reduces to a
nearest-neighbour interpolation scheme whose sole function is to send the irregularly positioned measure-ments to a uniform grid. The resulting reconstruction error scales with the distance between grid point and nearest data point, multiplied by the local gradient. It can be estimated using a known reference surface similar to the surface of interest, by comparing the known surface with its reconstruction from irregularly sampled data. Assuming that no other source of error is involved (completeness), and that these errors act independently (orthogonality), an overall root-mean-squared error can finally be estimated by summing the noise, attenuation and reconstruction errors squared, then taking the square root.
4.4 Validation tests
To check both the stereo algorithms and error estima-tion procedures, tests with a surface of known shape were conducted. The shape chosen is a wood model of a rippled bed, machined to high accuracy and rough-ened with fine sand, having a surface of sinusoidal elevation (see Table1 for precise parameters). A single camera combined with a mirror is used, in the config-uration shown in Fig.6. Fluorescent particles illumi-nated by black light dispersed onto the surface are used as markers for the stereo procedure (see Fig.7a). A series of 27 snapshots were acquired, with particles manually redistributed on the surface before each snapshot.
Maps of the true and reconstructed surfaces are presented in Fig.7. Panel 7b shows the sinusoidal sur-face of known amplitude and wavelength taken as ground truth. Panel 7c shows the surface obtained by sampling the true surface at the positions (xk, yk) of the
marker particles. This is what the reconstructed surface would look like if the stereo measurements of the par-ticle heights were perfectly accurate. Panel 7d shows the surface reconstructed from the actual stereo measure-ments using the maximum-likelihood approach of Sect.
4.2combined with the error balancing procedure of Sect.
4.3.
The different errors are plotted as a function of the binning subsample size n in Fig. 8. As explained in Sect.
4.3, three contributions can be distinguished, each characterised by a different dependence with respect to subsample size n: (1) noise due to random errors on the elevation of individual data points can be gradually averaged out by increasing n; (2) attenuation error, which grows as n increases; (3) these come on top of a reconstruction error, independent of n and estimated as the rms discrepancy between the true and reconstructed surfaces shown on panels b and c of Fig. 7. Note that errors (1) and (2) are estimated purely on the basis of the measured data, without using comparisons with the true surface. The thick black curve in Fig.8 represents the predicted overall error. On the same figure, crosses de-note the true error, evaluated as the rms discrepancy between the true and measured surfaces for various values of n. The value of n obtained by balancing noise and attenuation errors is n=12, and panel 7d shows the corresponding surface.
The predicted minimum rms error (i.e. the error ob-tained for the best choice of n) is 0.84 mm, while the
observed minimum rms error is 1.11 mm. The predicted error thus underestimates somewhat the error level. Various factors may account for this moderate offset: each error contribution is only approximately estimated; the different sources of error may interfere with each other; additional sources of error may play a role (e.g. departures from perfect perspective or calibration er-rors). Nevertheless, the predicted error level is close to the observed value, indicating that no significant error contribution was overlooked. The qualitative depen-dence of the error level on n is also well captured by the analysis. In particular, the predicted optimum sample size n is close to the observed optimum.
While not very precise, the proposed error estimation procedures are thus found to provide a valuable guide for the choice of bin size and for the estimation of the resulting error. This is especially important for the present experiments because, ultimately, measurements from two highly heterogeneous methods (stereoscopic and velocimetric) are to be compared with each other. As we experienced in the first stages of the present work, the application of different and arbitrary standards in designing the binning procedures for the two methods can lead to over-damping of one set of results and in-duce significant discrepancies between the two.
5 Velocimetric method
In contrast with the relatively direct character of the stereoscopic method, the velocimetric technique is indi-rect in nature. The technique seeks to reconstruct the surface topography on the basis of measurements of the horizontal velocity field. For that purpose, the method
Fig. 6 Imaging configuration for the stereo validation tests
exploits the Bernoulli relation applicable to free surface streamlines. This relation takes a particularly simple form if the flow can be approximated as a superposition of small perturbations upon a rapid mean current. To carry out the procedure, images from the top camera are first used to track particles on the water free surface and estimate their horizontal velocities. The elevation field is then obtained after suitable longitudinal and transverse
averaging of the individual velocity vectors. These steps are detailed in the next paragraphs.
5.1 Elevation-from-velocity principle
The following two assumptions are taken to apply lo-cally (on the scale of a few wavelengths) to antidune flows: (1) the time evolution of the loose sediment bed is sufficiently slow that it appears stationary to the rapidly flowing fluid, hence the flow can be taken as quasi-steady; (2) the flow can be considered inviscid (but it does not have to be assumed irrotational). The dynamics of the free surface can then be described by requiring conservation of mechanical energy along streamlines (see e.g. Batchelor 1967). Each surface streamline is subject to the following set of ordinary differential equations: dx dn¼ u ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi u2þ v2 p dy dn¼ v ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi u2þ v2 p dg dn¼ w ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi u2þ v2 p ; ð15a; b; cÞ d dn 1 2ðu 2þ v2þ w2Þ þ gg ¼ 0; ð16Þ
where n is a curvilinear coordinate measured along the streamline projection in the horizontal plane (x, y), u = (u, v, w) is the water velocity, and g is the free surface elevation. Equations 15a,b and c describe the streamline as a parametric curve [ x(n),y(n),g(n) ], along which the Bernoulli equation (16) is written. Vertical velocity w can be eliminated from the Bernoulli equation
Fig. 7 Stereo validation tests (run 0): a long exposure image of the dispersed particles, draped onto the test shape; b elevation map of the true sinusoidal surface; c reconstructed elevation map obtained by sampling the known surface at the particle positions; d maximum-likelihood elevation map reconstructed from the actual stereo measurements
Fig. 8 Error dependence on binning subsample size n for the stereo validation tests (run 0): predicted noise contribution (thin line), attenuation error (long dashes) and reconstruction error (short dashes); comparison of predicted (thick line) and measured (crosses) overall rms error
using (15c). It follows that if the horizontal components of the surface velocity field (u, v) of a steady flow are known (and if dissipation can be neglected), one can solve for the free surface elevation g along each streamline, up to a constant of integration. The constant can be determined by invoking a suitable boundary condition (point of known elevation along the stream-line) or stationarity condition (known average eleva-tion). This constitutes the basis of the velocimetric technique proposed here for the measurement of free surface topography. On a streamline-by-streamline ba-sis, the above principle can be applied to a wide class of flows beyond the present antidune case: examples in-clude flows over oscillatory bottoms, uniform currents past isolated accidents of topography, or shooting flows controlled upstream by a vane of known opening.
Before proceeding, it is useful to reduce the system to a dimensionless form. Given a characteristic scale [ u] for the velocity, a characteristic wave number can be defined as (Reynolds1965)
½k ¼ g
½u2: ð17Þ
Conversion to dimensionless variables can then be achieved by the transformation:
u! ½uu v! ½uv w! ½uw; ð18a; b; cÞ
x! x=½k y ! y=½k g! g=½k n! n=½k; ð19a; b; c; dÞ resulting in the following set of dimensionless equations: dx dn¼ u ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi u2þ v2 p dy dn¼ v ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi u2þ v2 p dg dn¼ w ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi u2þ v2 p ; ð20a; b; cÞ d dn 12ðu 2þ v2þ w2Þ þ g ¼ 0: ð21Þ
5.2 Approximation for small amplitude oscillations over a rapid mean current
An especially simple situation arises if the free surface can be approximated as small amplitude oscillatory perturbations superposed upon a rapid mean flow. This is the case for the present antidune flows, for which one can further assume that: (1) the mean flow is oriented in the longitudinal x direction; (2) the mean flow is approximately uniform, with only a weak mean velocity gradient along the transverse y direction. In that case, it is possible to decompose velocities into mean and per-turbation components and introduce a small parameter such that
uðx; yÞ ¼ u0þ u0¼ u0ðyÞ þ eu1ðx; yÞ
vðx; yÞ ¼ v0¼ ev1ðx; yÞ; ð22a; bÞ
where u0 is the mean current velocity, and e=[u¢]/[u]
measures the strength of the perturbation velocities u’, v’ (having characteristic scale [ u’]) relative to the mean current velocity u0(having characteristic scale [ u]). One
can then substitute expansions
g¼ g0þ e g1þ Oðe2Þ ð23Þ d dn¼ d dn 0 þe d dn 1 þOðe2Þ ð24Þ
in Eqs. 20 and 21 to obtain, at order 0,
g¼ g0; ð25Þ
and at order 1,
g1¼ u0u1 ð26Þ
where it was assumed that transverse variations in mean velocity ¶u0 /¶y are at most O(). One can finally
translate back to dimensional variables and obtain the key result g¼ g0þ g0¼ g0u0u 0 g þ ½u2 g Oðe 2Þ; ð27Þ
which is particularly simple both to use and to interpret: superimposed upon a mean surface level g0, the free
surface topography perturbations g¢ are proportional to the perturbations in longitudinal velocity u’, but oppo-site in phase. For torrential flows, the water velocity is slowest at the crests and fastest at the troughs, an observation pointed out already by Leonardo Da Vinci in his notebooks several centuries ago (Da Vinci ca.
1500/1975).
In practice, the ratio of fluctuation velocities to mean velocitiese=[u¢]/[u] obtained in the present experiments is about 0.1, with mean velocities of about 1 m/s. This leads to errors on the perturbed elevation of the order of (1 m/s)2(0.1)2/(10 m/s2)1 mm. More precise estimates obtained by calculating second order corrections are listed in Table1. These second order corrections are comparable in amplitude to measurement errors. Be-cause they are also highly sensitive to measurement noise, only the first order estimate is retained in what follows. Whatever the order of approximation, an obvious limitation of the velocimetric method is that only the topography perturbations are accessible, not the mean elevation. If needed, this has to be obtained by other means.
5.3 Horizontal velocities in world coordinates
Horizontal velocities (u,v) are obtained as follows from the image sequences acquired with the top camera (viewpoint V=T). First, particle positions on the digital images {RðTÞj;m } are pinpointed using the algorithm of Sect. 3.1. Particles are then tracked in the image plane using the Voronoı¨ match algorithm of Sect. 3.2. Index
pairs j(i)) are obtained, identifying matches {RðTÞj;m }, {RðTÞjð1Þ;mþ1 } likely to correspond to the same physical particles on successive images.
To convert image coordinates {RðTÞj;m } into world coordinates {ri,m}, the image positions are first projected
onto a plane corresponding to the mean water surface g0
(x,y). This is carried out as follows. To each particle image R(T) viewed under camera viewpoint T is associ-ated a ray
rðTÞ¼ pðTÞþ a qðTÞ; ð28Þ
where p(T) and q(T) are given by Eq. 3a and b. On the other hand, plane P 0can be specified by its Cartesian
equation a0 x+b0y+c0z=d0, where n0=(a0, b0, c0) is
the unit normal to the plane and d0=r0Æn0 its signed
distance from the origin (with r0an arbitary an arbitrary
point belonging to the plane). Upon substitution, the parameter value a corresponding to the point of inter-section of the ray with the plane is given by
a¼d0 p ðTÞ n 0 qðTÞ n 0 : ð29Þ
In the experimental configuration, plane P 0is close
to horizontal with Cartesian equation z=g0. Once world
coordinates ri,m=(xi,m, yi,m, zi,m) have been
approxi-mated by their projections in plane P 0, the horizontal
components (ui,m, vi,m) of velocity vector ui,m are
esti-mated from: ui;m¼ xjðiÞ;mþ1 xi;m Dt vi;m¼ yjðiÞ;mþ1 yi;m Dt ; ð30Þ
where Dt is the time interval separating successive frames.
The above procedure would be sufficient to obtain accurate measurements of the horizontal velocities (u,v) if either the water free surface was perfectly planar, or if we had been able to place the camera far above the water surface with a precisely vertical optical axis. In practice however, the free surface is perturbed by the antidune oscillations and the top camera is oriented slightly ob-liquely (as shown in Fig.2). As a result, a non-negligible bias results when the projection is performed onto a hypothetical non-perturbed plane. To correct for this bias, an iterative procedure is used: horizontal velocities derived from particle motions projected onto the mean water surface are used to construct a first estimate of the perturbed free surface topography. Particle motions are projected again onto the estimated perturbed topogra-phy, then a new horizontal velocity field is obtained, and so on. Since free surface slopes depart only moderately from horizontal, this procedure converges very quickly and two iterations were found sufficient.
5.4 Averaging of the velocity field
An averaged horizontal velocity field remains to be ex-tracted from the full set of horizontal particle velocities
(ui,m, vi,m) defined at horizontal positions (xi,m, yi,m).
Rather than combining right away the vectors corre-sponding to different times tminto one single set to be
indiscriminately binned and averaged, the regularisation procedure adopted here seeks to preserve the streamline structure which underlies the various individual velocity vectors. This structure presents important implication for the error formation process, and it is sought to ac-count for it by splitting the averaging procedure into two steps: (1) along-trajectory averaging; (2) transverse averaging. For each of the two steps, a balance between attenuation and noise error can again be sought in order to minimise the combined error.
Introducing hats to distinguish measured estimates from hatless ‘‘true’’ quantities, consider measurements ^
um acquired at successive time instants tm along a
single particle trajectory. We have the standard par-ticle tracking velocimetry (PTV) estimate (e.g. Adrian
1991) ^ um¼ ^ xmþ1 ^xm Dt ; ð31Þ
where measurement ^xm approximates the true particle
position x(tm) at time tm, and measurement ^um
approx-imates the true instantaneous velocity u(tm+1/2) at
intermediate time tm+1/2=tm+(1/2) Dt. A simple
sto-chastic model for error formation can be derived by assuming (Wernet and Pline 1993) that measured posi-tions ^xmcorrespond to the actual particle position x(tm)
to which is added Gaussian white noise, i.e.
^
xm¼ xðtmÞ þ rxcm; ð32Þ
where rxis the root mean square (rms) error on particle
position, and where the cm are independent Gaussian
random variables of zero mean and unit standard devi-ation. Invoking (31) and (32) as well as Lagrangian relation x(t) = u(t)d t, velocity estimate ^um can be
expressed in terms of the actual velocity history ut) as ^ um¼ umþ rx Dt cmþ1 cm ; ð33Þ where um¼ 1 Dt Z tmþ1 tm uðtÞdt: ð34Þ
This constitutes an observation equation (Honerk-amp 1998) for the PTV measurement process, i.e. an idealised description of the connection between the true signal and the measured one. From expression (33), two discrepancies are seen to arise between measurements ^um
and the actual instantaneous velocities ut). The first discrepancy is a corrupting high frequency noise due to random position error. The second discrepancy results from averaging over interval Dt, which attenuates high frequency velocity variations.
Minimisation of the combined error can again be sought by balancing the effects of noise and attenuation.
A filtered signal ~uðnÞm can be obtained by convoluting signal ^umwith a rectangular window of width (2 n+1)Dt:
~ uðnÞm ¼ 1 2nþ 1 Xn l¼n ^ umþl: ð35Þ
(The same operation can also be applied to the transverse velocity measurements ^vn to yield
along-tra-jectory filtered values ~vðnÞm :) Substitution of (35) into (33) yields an observation equation for the filtered velocities
~ uðnÞm ¼ uðnÞm þ rx ð2n þ 1ÞDtðcmþ1þn cmnÞ; ð36Þ where uðnÞm ¼ 1 ð2n þ 1ÞDt Z tmþ1þn tmn uðtÞdt: ð37Þ
Relation (36) is equivalent to the earlier observation Eq. 33, except that it is applied to the longer time interval (2 n+1)Dt. From Eq. 36, the rms error on the filtered velocity due to position error is given by
ð~uðnÞm uðnÞm Þ2 D E1=2 ¼ ffiffiffi 2 p rx ð2n þ 1ÞDt: ð38Þ
Because measurements are obtained along a trajec-tory, position errors cancel each other along the way and only the endpoint errors matter. As a result, the rms error on velocity (Eq. 38) due to position errors creases in proportion to 1/(2 n+1), a much faster de-crease with filter width than the factor 1=pffiffiffiffiffiffiffiffiffiffiffiffiffiffi2nþ 1;which would be obtained if errors on velocity were uncorre-lated Gaussian random variables. It is to account for this feature of the error formation process that a two step procedure starting with along-trajectory filtering is selected in preference to indiscriminate binning of the whole set of velocity measurements. Because the present experiments feature variations of relatively long wave-length, one can approximate uð1Þm uð2Þm ;yielding for the random position error rx the rough estimate
rx¼
3 5 ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2ð32þ 52Þ
p Dð~uð1Þm ~uð2Þm Þ2E1=2: ð39Þ The resulting values for the different runs are listed in Table1. Converted back to pixel coordinates, the par-ticle positioning error is around 0.1 pixel. Due to the limited resolution of the camera, such small errors in pixel units still translate into significant errors on velocities.
The corresponding attenuation error can be esti-mated by idealising the velocity history of a particle along its streamline as a harmonic signal, i.e.
uðtÞ ¼ u0þ a sinðx tÞ; ð40Þ
where angular frequencyx 2p[ u]/[k], amplitude a [u¢], and characteristic scales [ u], [k] and [ u’] can be estimated for a given run (see Table 1). Alternatively,
Eq. 17 can be invoked to approximate [k ]=2p /[k] 2p [u]2 /g. Convoluting signal (Eq. 40) with a rectangular window of width (2 n+1)Dt results in filtered signal
uðnÞðtÞ ¼ u0þ aðnÞsinðx tÞ; ð41Þ
where the damped amplitude aðnÞis given by (e.g. Ja¨hne
1995) aðnÞ
a ¼ sinc(12xð2n þ 1ÞDtÞ ð42Þ
in which again sinc(x) = sin(x)/ x. Developing Eq. 42 by Taylor around the origin and setting a=[u¢], the rms attenuation error on the filtered signal can be approxi-mated by ðuðnÞ uÞ2 D E1=2 ¼ 1pffiffiffi2ða aðnÞÞ p 2 6pffiffiffi2 ½u0½u2ðDtÞ2 ½k2 ð2n þ 1Þ 2 : ð43Þ
Balancing position noise error (38) and attenuation error (43), one obtains for the conditions of the present experiments an optimal filter width (2n+1)=13. The estimated error on trajectory-filtered velocities due to either noise or attenuation is of the order of 0.01 m/s.
The last remaining operation is transverse averaging. This is performed in two steps. First, filtered along-tra-jectory velocities are sampled at the intersections of the trajectories with transverse planes placed at equal intervals xialong the x-axis. This yields for each xia set
of irregularly distributed data points ðyk; ~uk; ~vkÞ: The
second step of the procedure consists of converting these measurements into regularly spaced data at positions yj
of a uniform grid. The averaging is performed by way of overlapping bins of equal width Dy, where the bin width Dy is again chosen to minimise the combined error resulting from noise and attenuation. Noise at this stage comes both from the imperfect PTV measurements and from fluctuations in time of the physical trajectories, and is assumed to be Gaussian. Using the same error bal-ancing approach as before (see Sect.4.3), the optimal bin width is found to be Dy35 mm. The estimated error due to either noise or attenuation at the transverse averaging stage is again of the order of 0.01 m/s. Noise and attenuation errors incurred at both the streamline and transverse averaging stages yield an overall rms velocity error of the order of 0.02 m/s. This amounts to errors of the order of 2% relative to the mean velocity (see Table1). Upon invoking the approximate Bernoulli equation (27), elevation errors due to velocity measure-ment errors are of the order of 2 mm.
To these must be added the neglected contribution from second order perturbations. Assuming the errors to be independent, an overall error is again estimated by adding the squares of the rms contributions and taking the square root. The resulting values for the various runs are given in Table1. Reconstruction errors arising from
interpolation onto a regular grid are not significant for the velocimetric method because interpolation is con-ducted to first order accuracy. The interpolants used are piecewise linear rather than piecewise constant. The latter approach, less accurate but more robust, is used for binning in the stereo method due to the presence of stereo mismatches.
6 Results and discussion
Examples of free surface topography measurements obtained using both the stereoscopic and velocimetric methods are presented in Figs.9,10,11, 12, 13and 14. Taken from a series of ten runs, the three experimental runs shown (runs 3, 5 and 6) exhibit a range of different patterns: wide-crested rolls for run 3 (Fig.9), a zig-zag pattern of staggered peaks and troughs for run 5 (Fig.10), and a narrow wave train with sharply-peaked crests on one side of the channel for run 6 (Fig.11). The gray-coded surface elevation maps of Figs.9,10and11
provide full-field pictures of these different patterns. Since the velocimetric method cannot determine the mean water level g0, the latter has been subtracted from
the stereo results to obtain for both methods the per-turbed topography g¢(x,y)=g(x,y) g0.
These maps allow qualitative comparisons between the two methods. Some superficial differences between the stereoscopic and velocimetric results can be noted. First, due to the different binning procedures used, the spatial coverage of the two methods does not fully coincide. Because the corresponding binning assigns an elevation even to areas where few particles were identi-fied, the stereo results extend over the whole domain. The velocimetric results, on the other hand, cover only the restricted domain where particle trajectories were retrieved. Second, the surfaces exhibit distinct textures that reflect the peculiarities of each of the two methods. Because the density of the cloud of three-dimensional points can vary abruptly, the maximum-likelihood fit-ting used in the stereo method leads to staircase-like surfaces. On the other hand, the velocimetric maps are smoother, but feature longitudinal stripes due to the streamlines along which velocities are measured and averaged. These are fine-grained effects, however, which do not affect the overall surface elevation maps.
Qualitatively, the agreement between the two meth-ods is seen to be quite good for all three of the patterns examined. The magnitudes, locations and arrangements of the peaks and troughs are in good correspondence.
Fig. 9 Reconstructed free surface topography g¢(x,y) for run 3 (broad crested rolls): a stereoscopic method; b velocimetric method
Fig. 10 Reconstructed free surface topography g¢(x,y) for run 5 (zig-zag pattern): a stereoscopic method; b velocimetric method
Local regions of discrepancy are nonetheless present, for instance the spurious crater-like feature exhibited by the stereo results of Fig. 11a near position (x,y)=(150,400), but absent from the velocimetric results of Fig.11b. This is a region of wave breaking where the stereo method
appears to have failed, likely due to an insufficient number of successfully matched stereo pairs in that zone.
To facilitate precise comparisons, the same results are plotted again in a different format on Figs.12, 13 and
14. Here profiles g¢(x) reconstructed by the two methods are plotted together for various positions y across the width of the channel. The stereo profiles (thin lines) and velocimetric profiles (thick lines) are found to be in reasonably close qualitative and quantitative agreement. Results in the central region of the viewing volume are generally better than the results near the edges. Espe-cially poor results are recorded for the stereo profile of Fig.14a, which passes through the ‘‘crater’’ zone of Fig.11a. For this run (run 6), the stereo cameras were placed lower with respect to the surface than for runs 3
Fig. 14 Comparison of free surface profiles for run 6 obtained by the stereoscopic (thin lines) and velocimetric (thick lines) methods at selected sections: a y=400 mm; b y=300 mm; c y=200 mm; d y=100 mm
Fig. 11 Reconstructed free surface topography g¢(x,y) for run 6 (narrow wave train on one side of the channel): a stereoscopic method; b velocimetric method. Circled: crater-like region where stereo measurements break down
Fig. 12 Comparison of free surface profiles for run 3 obtained by the stereoscopic (thin lines) and velocimetric (thick lines) methods at selected sections: a y=400 mm; b y=300 mm; c y=200 mm; d y=100 mm
Fig. 13 Comparison of free surface profiles for run 5 obtained by the stereoscopic (thin lines) and velocimetric (thick lines) methods at selected sections: a y=400 mm; b y=300 mm; c y=200 mm; d y=100 mm
and 5, leading to near-occlusion in the troughs in addi-tion to the difficulties associated with wave breaking at the peaks.
Computation of the root-mean-squared discrepancy Æ(g¢stereo g¢vel)2 æ1/2 between the results of the two
methods (averaged over a square domain of 500·500 mm) for the three runs 3, 5 and 6 yields rms errors of 3.8, 5.1 and 4.7 mm, respectively. Relative to elevation range g¢max-g’min40 mm, this amounts to
relative errors of the order of 10–15%, which is com-parable to the error level encountered in the stereo val-idation tests of Sect.4.4(see also Table1). If instead one gauges errors against the wavelength, taken as repre-sentative macroscopic length scale, relative errors of the order of 1–1.5% are obtained. The following rule of thumb thus appears applicable here: to obtain 10% accuracy for the perturbed velocities and elevations, 1% accuracy relative to the mean flow quantities (mean velocity and wavelength) is required.
The detailed values of the predicted error contri-butions are listed in Table 1 for the three different runs. To estimate the stereo reconstruction error, the surface obtained by the velocimetric technique was used as reference surface. The observed discrepancies (of the order of 4–5 mm) turn out to be somewhat higher than the predicted discrepancies (of the order of 2–3 mm), obtained by summing the squares of the predicted error levels in both techniques, then taking the square root. This slight underestimation of the error levels may be ascribed to three factors: (1) the approximate nature of the error analysis procedures; (2) model errors (for instance deviations of the image formation process from a perfect perspective projec-tion, or departures of actual streamline dynamics from the assumed Bernoulli relation); (3) the presence of regions where one of the two methods performs par-ticularly poorly, as for instance in the ‘‘crater zone’’ of Fig.11a mentioned previously.
Overall, the quality of the comparison is encouraging. Errors are not negligible, but remain within reasonable bounds considering that both methods yield whole-field measurements rather than point values. Furthermore, the order of magnitude of the incurred errors appears consistent with the error estimates derived using the techniques of Sects. 4and 5. Errors not included in the analysis can therefore be ascertained to be small. Despite incipient wave breaking, for instance, streamline dynamics appear well approximated by the Bernoulli relation on the scale of a few wavelengths. Note that this implies only that the flow is nearly inviscid, since the Bernoulli equation applied on a streamline-by-stream-line basis does not require the flow to be irrotational (see e.g. Batchelor 1967).
Most importantly, the spatial patterns associated with flows over antidunes are successfully captured by both methods. The results vividly depict a variety of motifs: broad crested rolls, zigzags and narrow train of peaks and troughs. Due to the complex way in which light interacts with the rough water surface, these
organised patterns could not be so easily grasped by pure visual inspection of the laboratory flows.
7 Conclusions
In the present work, two approaches to the measurement of free surface topography in high Froude number flows were presented and compared. Both rely on the imaging of floating tracers dispersed on the surface, but they resort to different reconstruction principles: stereoscopic matching of rays associated with two camera viewpoints, on the one hand, and indirect estimation of the eleva-tions based on the horizontal velocity field measured using a single camera, on the other hand. A waveform of known shape was used to check the stereo algorithms and error estimation approach. The results of the two techniques were then compared for flows over antidunes featuring various surface patterns. Although somewhat larger than expected based on error estimates, the ob-served discrepancies were found to be reasonably small. The obtained relief maps vividly depict the variety of motifs that can evolve as a result of interaction between shallow flows and loose sediment beds.
Drawbacks, advantages and possible improvements can be identified for both techniques. Stereoscopic matching and velocimetric tracking both involve pairing particles on distinct images, but the stereo methods were found to generate a much greater proportion of mis-matches (up to 40% of the data points). This stems partly from the fact that the stereo technique matches isolated particles whereas the velocimetric technique tracks patterns of neighbouring particles. The higher vulnerability of the stereo technique to mismatches has various consequences: greater accuracy requirements for the camera calibration procedure; tighter limits on the density of particles which can be reliably imaged; pre-cedence of robustness over accuracy when interpolating data. Camera synchronisation is also a paramount concern, best addressed by synchronising two sensors at the time of image acquisition (rather than a posteriori as in the present work).
A key advantage of the stereo technique is that, in principle, it applies to general flows independently of any specific assumption (occlusion being the only limi-tation). In practice, the stereo technique requires the flow to vary slowly in order to combine data from multiple image pairs and obtain a surface of sufficient resolution. The velocimetric technique, by contrast, fundamentally depends on certain physical hypotheses. The assumed Bernoulli equation requires the flow to be quasi-steady, and its simplified version further assumes moderate perturbations superimposed upon a quasi-uniform mean current. It also requires the mean surface level to be known by other means. When these condi-tions are met, however, the present results demonstrate that the velocity-based approach works well, producing reasonable topography measurements. For the present experiments, estimated error levels for the velocimetric
