Stereoscopic and velocimetric reconstructions of the free surface topography of antidune flows

4. Stereoscopic reconstruction

The first reconstruction method is a feature-based stereoscopic matching procedure used in digital photogrammetry (see e.g. Wolf and Dewitt 2000), 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 developed 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 and on images seen by the left and right cameras (V = L, R). The corresponding rays are given by parametric equations (see eq. 2)

, .

(4a,b)

With reference to Fig. 4, the points on the left and right rays corresponding to the minimum distance between the two rays are parameterised by the values , which are solutions to the system

,

(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 midpoint and distance between the two points of closest encounter are then given by

; .

(6a,b)

Figure 4. Intersection of two rays associated with the left and right camera viewpoints.

Due to the limited accuracy of the camera calibration and image plane measurements, rays corresponding to the same physical point will not perfectly intersect. Midpoint constitutes an approximate intersection, while the distance of closest encounter provides an indicator of the quality of the approximation.

A procedure to select stereo pairs follows. Let and designate the two sets of particle positions picked up at the same time t  by the left and right cameras. The task is to find the pairing j(i) 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 sketched earlier in section 3.2 and 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

(7)

under the constraint that a given particle image can participate in only one stereo pair. This bijective 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 . 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 acquisition, 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 section 3.2. Interpolation in time is then used to resample their positions at the instants 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 cameras. This problem could be solved by using a stereo set-up 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 {r_i,m} positioned within the viewing volume at time instants {t_m}. Since in the present application the flow is assumed to be quasi-steady over the duration of each video acquisition, the different subsets {r_i,m} can further be pooled into one single set {r_k}. A continuous water free surface elevation η(x, y) remains to be extracted from this cloud of points. If the matching process were perfectly 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 errors. The measured point set {r_k} includes two distinct subsets: (i) 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; (ii) points resulting from mismatched rays. It is reasonable 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.

Figure 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).

With reference to Fig. 5a let us suppose for the sake of argument an idealised viewing geometry with focal points of the projection sent to infinity, and a uniform distribution of particles on a wavy surface. Geometrically, mismatches correspond to the establishment of a wrong correspondence between two rays that approximately 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 approximate 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 Δx, Δy around position (x, y). Let ω ₁ designate the correctly matched state in which a particle can find itself, and ω ₂ designate the mismatched state. P(ω ₁) and P(ω ₂) are the corresponding probabilities. Based on the above model, the point elevations {z_k} are randomly distributed with probability density function given by

(8)

where Θ₁, Θ₂ designate the parameter vectors which characterise the distribution of each class. This is a so-called mixture density, in which conditional densities , are called component densities and probabilities , are called the mixing parameters (see e.g. Duda and Hart 1973). For the normal and uniform distributions postulated above, we have

,

(9a,b)

hence parameter vectors are given by Θ₁ = (η, σ_η ) and Θ₂ = (z_a, z_b), where η, σ_η  are the mean and standard deviations of the distribution of correct matches and z_a, z_b are the (possibly truncated) vertical limits of the viewing volume at location (x, y). Symbol η has been used to designate the mean of the correct matches because 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 η 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 z_a and z_b are known and, by definition, . Three parameters thus remain to be determined: η, σ_η  and . In a manner analogous to the case of normal mixtures treated by Duda and Hart (1973), maximum-likelihood estimates , and must satisfy

,

(10a)

,

(10b)

,

(10c)

where, by Bayes' rule (see Duda and Hart 1973)

(11)

and where and are given by distributions (9). Despite the involved expressions, estimates (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 , and can first be used to evaluate 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 guaranteed in all cases, the procedure was found to perform adequately in this application, provided reasonable initial estimates are chosen. For the present antidune experiments, the converged estimates for and are  ≈ 3 mm and  ≈ 0.6. The latter number is especially significant: about 40 percent 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.

Figure 5b 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 distributed 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 clustered points. Similar techniques aimed at weighting down outlying measurements in the presence of non-Gaussian noise are used in high-energy physics (e.g. Frü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 x_i× y_j. 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 following considerations. When applied to a spatially distributed set of data corrupted by uncorrelated Gaussian noise, averaging produces two different effects: on the one hand, the Gaussian noise error diminishes in proportion to 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 attenuation 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 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

,

(12)

where the triangular brackets denote an ensemble average, is the averaged surface which would be obtained by bin-averaging in the absence of noise, is the standard deviation of the Gaussian error as estimated from eq. (10c), and is the number of Gaussian distributed points estimated from eq. (10a). Reduced estimate 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 .

An estimate of the attenuation error, on the other hand, can be obtained as follows. Let ν 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 number of points n to each bin is given by . Assume now that the free surface is roughly harmonic, with variations in the x and y directions characterised by wavelengths λ_x ≈  λ_y ≈ [ λ] and amplitude a ≈ [η'], where [λ] and [η'] designate characteristic horizontal and vertical scales of the free surface oscillations. The averaging process can be idealised as a convolution of the harmonic surface with a square top-hat filter of size Δx × Δy. This results in an attenuated amplitude given by (e.g. Jähne 1995)

,

(13)

where the sinc function is defined as sinc(x) = sin(x)/x. Developing the sinc function by Taylor around the origin (up to second order), one can estimate the attenuation error due to the spatial averaging as

,

(14)

where the attenuation error is seen to grow in proportion 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 Table 1), this results in an optimal number of points n ≈ 50 and an average bin size Δx = Δy ≈ 45 mm. The associated rms error on surface elevation due to either noise or attenuation is around 0.5 mm.

Table 1. Experimental parameters and error estimates.

In addition to these two sources of error which depend 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 (), the procedure reduces to a nearest-neighbour interpolation scheme whose sole function is to send the irregularly positioned measurements 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 estimation 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 roughened with fine sand, having a surface of sinusoidal elevation (see Table 1 for precise parameters). A single camera combined with a mirror is used, in the configuration shown on Fig. 6. Fluorescent particles illuminated 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.

Figure 6. Imaging configuration for the stereo validation tests.

Maps of the true and reconstructed surfaces are presented on Fig. 7. Panel 7b shows the sinusoidal surface of known amplitude and wavelength taken as ground truth. Panel 7c shows the surface obtained by sampling the true surface at the positions of the marker particles. This is what the reconstructed surface would look like if the stereo measurements of the particle heights were perfectly accurate. Panel 7d shows the surface reconstructed from the actual stereo measurements using the maximum-likelihood approach of section 4.2 combined with the error balancing procedure of section 4.3.

Figure 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.

The different errors are plotted as a function of the binning subsample size n on Fig. 8. As explained in section 4.3, three contributions can be distinguished, each characterised by a different dependence with respect to subsample size n: i) noise due to random errors on the elevation of individual data points can be gradually averaged out by increasing n; ii) attenuation error, which grows as n increases; iii) 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 i) and ii) are estimated purely on the basis of the measured data, without using comparisons with the true surface. The thick black curve on Fig. 8 represents the predicted overall error. On the same figure, crosses denote 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 , and panel 7d shows the corresponding surface.

Figure 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.

The predicted minimum rms error (i.e. the error obtained 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 errors). Nevertheless, the predicted error level is close to the observed value, indicating that no significant error contribution was overlooked. The qualitative dependence 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 induce significant discrepancies between the two.