High-resolution single-camera photogrammetry: incorporation of refraction at a fluid interface

Photogrammetry uses images of a three-dimensional structure to derive information about its shape and position. In this work, a photogrammetric technique is implemented with a single camera and a digital projector to measure changes in an underwater sediment bed. This implementation incorporates refraction at an interface allowing for measurements through a deformed or changing water surface. The digital projector provides flexibility in choosing projected patterns and has a high frame rate, which allows to easily increase the spatial and temporal resolution of the measurements. The technique requires first for both the camera and the projector to be calibrated using triangulation. With the calibration, we construct lines in three-dimensional space that originate from the projector and the camera, and intersect on the surface to be measured. To correctly incorporate refraction due to a change of medium, each line in space is recalculated from its intersection with the interface using Snell’s law. This has the benefit that only one calibration for measurements is needed if the location and shape of the interface are known. The technique is validated by measuring a submerged undulated surface, plastic objects and a sediment bed. In particular, the undulated plate is reconstructed under a flat and a parabolic water surface. Finally, the technique is used in combination with particle image velocimetry to dynamically measure a changing sediment bed under an oscillating flow and the flow velocity at the free surface.


Introduction
The development and continued advancement of digital technologies have led to the improvement of digital imaging, increasing image resolution, acquisition speed, memory storage and decreasing transfer time. This has helped 3 Page 2 of 19 advance a variety of research areas and applications that use two-dimensional images as a foundation for both qualitative and quantitative analyses. One such application is photogrammetry, which uses multiple images of a subject (each from a different position) to derive information about its three-dimensional structure. This is done using perspective projection (Longuet-Higgins 1986) or triangulation (Sansoni et al. 1997), where knowledge of the orientation of the image planes and the position of the imaging device (camera) are used to determine the position and elevation of the object (or surface) in the image.
Photogrammetry has been used in a wide variety of applications such as topographic modeling (Remondino and El-Hakim 2006), morphometric measurement of corals (Bythell et al. 2001), archeological surveys (Henderson et al. 2013;McCarthy and Benjamin 2014), and several topics in geomorphology (Haneberg 2008;Keutterling and Thomas 2006;Lane 2000;Marzolff and Poesen 2009;Pyle et al. 1997;Westaway et al. 2003). Furthermore, since photogrammetry is non-intrusive, it is capable of measuring very delicate structures at different scales, such as solar sails and antennas (Giersch 2001;Pappa et al. 2001) and flapping microwings (Curtis 2009). Within the study of fluid motion, photogrammetry has been used for surface flow visualization (Karthikeyan and Venkatakrishnan 2011) and for measurements of the shape of free surface water waves (Rupnik et al. 2015). The capability of photogrammetry to measure surfaces in a non-intrusive way makes it a good alternative to quantify changes in sediment beds in laboratory experiments (Baglio et al. 2001;Peter Heng et al. 2010;Rieke-Zapp and Nearing 2005;Stancanelli et al. 2011;Stojic et al. 1998).
Photogrammetry competes with several other techniques (each with its own advantages and limitations) due to the ample range of applications of sediment bed morphodynamics for coastal and river engineering, and in fundamental research. For example, despite the poor spatial resolution, point gauges (Lança et al. 2013) have been thoroughly used to measure the sediment bed thickness due to the high temporal resolution. The devices used for this type of measurements range from photosensors (Ballio and Radice 2003) and ultrasonic sensors (Best and Ashworth 1994;Coleman et al. 2003), to infrared beams (Richards and Robert 1986). Although the spatial resolution can be improved using multiple point gauges, it is often costly and intrusive (Sheppard and Miller Jr 2006). It is also possible to measure the sediment bed thickness using electrical field attenuation (De Rooij et al. 1999) or light attenuation (Munro and Dalziel 2005), but these methods require complex, specific measurement setups, for example, a particular type of particles. Compared with planar laser sheet illumination (Younkin and Hill 2009), photogrammetry has the clear advantage of being able to measure the bed thickness over an area and no real disadvantage in terms of accuracy ). Photogrammetry has been compared with two more traditional techniques (calibrated pile and echosounder) to measure the scour hole around a pile (Porter et al. 2014). In this case, photogrammetry proved to be superior in terms of accuracy to the echosounder for a much more affordable price. The advancement of digital imaging allows photogrammetry to provide measurements with increasingly higher spatial and temporal resolution, which is ideal for measuring evolving sediment beds.
The photogrammetric techniques used to measure changes in a sediment bed consist of a projection device and at least one acquisition device (a camera). The use of two cameras in a stereoscopic array is also possible Peter Heng et al. 2010;Sumer et al. 2012). The projection device, which is either a digital projector or a laser, projects a certain pattern (usually composed of dots) onto the bed. The use of light patterns reduces the errors related to the positioning of the measuring points in comparison to feature recognition . The changes in the projected pattern due to changes in the surface are registered by the camera. Changes in a free water surface have also been measured with this technique (Cang et al. 2019;Gomit et al. 2013).
Recent studies have used photogrammetry to measure changes in the morphology of a sediment bed inside a flume (Butler et al. 2002;Stancanelli et al. 2011). However, measuring through the water interface is a common problem, which depends on the amplitude of the disturbances and on the distance from the measurement plane to the interface (as discussed by Akutina et al. 2018, for particle tracking velocimetry measurements through a free surface ). Dynamic measurements were obtained by Baglio et al. (2001) but did not account for distortions produced by the lens or effects due to refraction. A similar system was implemented in a study by Sumer et al. (2012), where corrections for lens distortions were included but refraction was not modeled correctly. Analytical, geometry-based corrections have also been implemented to account for refraction by Butler et al. (2002) and Dietrich (2017). A solution to reduce the effect of refraction is to use underwater cameras Porter et al. 2014), but this is not always possible. In a more recent study, the refraction at the water interface was taken into account during the calibration while using a pattern composed of parallel lines to measure underwater objects and the thickness of a sediment bed (Cebrián-Robles and Ortega-Casanova 2016).
In the current paper, a version of the photogrammetric technique using a single camera and a single digital projector is proposed to measure the morphology of a sediment bed through the air-water interface. We exploit the advantages that using a digital projector gives: flexibility when choosing the patterns, the size and distribution of the measuring points, the number of patterns projected and their rate of projection. To incorporate the effects of refraction through an interface, we use Snell's law. This is done in such a way that a single (dry) calibration for both the camera and the projector can be used for measurements under any water height provided that the position of the interface is known. Furthermore, this calibration also allows to take other optical effects (such as lens deformations) into account.
This version of the measurement technique has then four main advantages: (1) it allows for measurements through an interface (e.g., a free surface), (2) the flexibility gained from using a digital project allows for high temporal and spatial resolutions, (3) to easily incorporate other measurements, such as flow velocity measurements through particle image velocimetry (PIV), and (4) to tailor the technique to a wide variety of experimental setups. However, the technique has certain limitations as well. For example, the position of the interface has to be accurately known, and measurements have to be performed under clear-water conditions.
To show the robustness of the method, surface measurements are performed under different test conditions. For this, two different experimental setups were used. In the first setup, complex objects and a stationary sediment bed are measured under a flat and a parabolic air-water interface. In the second experimental configuration, we measured the sediment bed height distribution under an oscillating flow and oscillating water surface height. For this latter case, the flow velocities are simultaneously measured with PIV.
The present paper is organized as follows. A description of the developed method is given in Sect. 2. In Sect. 3, the first experimental setup is described, and results of the reconstruction of the complex objects and surfaces are presented. The setup and the results for the oscillating flow are presented in Sect. 5. Results are discussed in Sect. 6, and the conclusions are outlined in Sect. 7.

Method
The measurement technique uses a digital projector set at an angle and a single overhead camera that looks vertically down. The projector is used to project a series of images onto the surface to be measured. The images consist of a number of elliptical bright patches over a black background. When the images are projected, each patch generates a light ray, originating from the projector and intersecting with the illuminated surface. The position of these intersections depends on both the angle and position of the projector and the position of the surface. The size of the elliptical patches and the projector's angle are such that these patches are viewed (and recorded) by the overhead camera as circular dots when projected onto a horizontal surface. Each dot has a pixel diameter and are separated by a distance e , as shown in the projected image Fig. 1a and on the recorded image by the camera Fig. 1b. Depending on the height distribution of the illuminated surface, the projected image will deform from the perspective of the camera.
The local surface height at a dot position can be determined by relating this dot in the recorded images to their respective originating light ray. This is equivalent to finding the intersection between the light ray going from the projector to the surface and the light ray going from the dot to the camera. A two-dimensional schematic of this is shown in Fig. 2. Since each dot is a point measurement, a three-dimensional surface reconstruction is made using all the measurements obtained from all ray intersections. In this way, increasing the number of dots per image (reducing e ) increases the resolution of the reconstruction, but this has limits since increasing the density of dots per image can lead to dot overlapping. Another option is to increase the number of projected images and combine their results. , projecting a set of dot patterns onto an arbitrary surface with a height h(x, y). The camera is placed at a height H C , looking vertically down. As an example, two projected dots are indicated by the black dots. The position of these projected dots without the surface h is indicated by the red dots. With the camera, photographs of the projected patterns are obtained. The blue lines indicate the light rays originating at the dot and captured by the camera, whereas the dotted lines indicate the light rays without the surface h. The black lines indicate the light rays from the projector Since the light rays obtained from both the camera and the projector are based on images, the location of the projected dots and measured height is given in pixels (X, Y). To transform this information into real-world coordinates (x, y), both the camera and the projections are calibrated. This also accounts for the deformation of the images caused by the projection angle, the focal length, and the warping produced by the camera lens.

Camera calibration
The camera calibration is done using images of a surface with known dimensions and coordinates, such as a square grid. In our case, a flat 60 cm × 60 cm plate (from here on referred to as the calibration plate) is photographed in the same horizontal position at different heights. The plate has on top an equidistant array of 23 × 23 white circular dots, each dot with a diameter of 0.5 cm and a 2.5 cm spacing between them. As shown in Fig. 3, even though the plate is only displaced vertically, the number of pixels between each white dot changes. This means that the relation between pixel and real-world coordinates is height dependent, such that x = F(X, Y, h) and y = G(X, Y, h) . Using a third-order polynomial fit (Soloff et al. 1997), we relate the pixel location (X, Y) of the white dots in the calibration plate at each height the plate was placed, to their real-world position (x, y). In other words, a surface fit is made for each height h i , where the subscript i indicates the different number of heights the plate was placed.
To determine the real-world location associated to a pixel for any height, we take a virtual square grid in pixel space, where the position of the grid points are (X l , Y l ) . The real-world coordinates of these pixel positions at every calibration height (x l , y l , h i ) are obtained using the previously calculated third-order polynomial fit. Afterwards, using the real-world coordinates of each point of the virtual grid at all calibration heights h i , lines are constructed using least squares. Each line passes through the real-world location of a virtual grid point and its corresponding locations at all different heights h i . These lines are analogous to the light rays from optical geometry. A two-dimensional schematic of this approach is shown in Fig. 4.
By extending the calculated lines in space to a predefined plane h = h 0 (typically the bottom of the domain so that h 0 = 0 mm ), the virtual grid is defined at the lowest plane ( X 0 , Y 0 ). Each line in space is then redefined as an euclidean vector �� ⃗ C = [dx c , dy c , dz c ] , with its origin at ( x c,0 , y c,0 , h 0 = 0 mm ). This vector is made non-dimensional and normalized in the vertical so that dz c i = 1 . This means that, for each pixel k, there is a vector �� ⃗ C k with its associated components and origin. Each of the components of �� ⃗ C k is determined using a third-order polynomial fit for X and Y and is strictly valid within the calibration area, but it can be extended outside the calibration area as long as the fit used remains within the accuracy set by the user. An example of a calibration is shown in Fig. 5, where the lines in space are obtained from the fit through eight planes. For a correct calibration, the lines in space intersect with the virtual grid points and converge in a small region where the camera is located.

Digital projector calibration
In a similar way to the camera calibration, the digital projector's calibration requires capturing a series of images projected on horizontal planes at different heights h i . At each plane height, all the images used to perform the measurement need to be projected and photographed. For given images, changes in the dot position are only due to the height difference. To determine the location of each dot, either an ellipse or a circle can be fitted to each bright patch to find a centroid, transforming each image to a list of (X, Y) positions. Although the shape of the dot depends on the height distribution of the surface onto which it is projected, we assume that the center of the dot remains relatively unaffected.
Using the camera calibration, the pixel location of all projected dots in an image ( X p , Y p ) is translated to realworld coordinates ( x p , y p ) for each height h i . Analogous to the procedure in the camera calibration, these real-world coordinates are used to construct lines in space using a linear regression method (least squares). These lines are then treated as vectors � ⃗ P = [dx p , dy p , dz p ] , with an origin in ( x p,0 , y p,0 , h 0 = 0 mm ). Here, ( x p,0 ,y p,0 ) is the horizontal Fig. 3 Photographs of the calibration plate at two different heights: a h = 1.5 cm and b h = 19.5 cm . The plate consists of a 23 × 23 grid of white dots, each with a diameter of 0.5 cm and separated by 2.5 cm from center to center real-world position of the projected dot at the bottom of the tank. Similar to the camera vector �� ⃗ C k , each projected dot j has a vector � ⃗ P j associated to it which is non-dimensional and normalized along the vertical, i.e. dz p = 1.
An example of lines in space constructed from the images composed of randomly placed dots projected at different heights is shown in Fig. 6a. For a correct calibration of the digital projector, the constructed lines in space pass through their corresponding projected dots and converge in a small region in space, which coincides with the position of the projector when modeled as a pinhole source. Figure 6b shows in a two-dimensional configuration the superimposed lines in space for both the camera and the digital projector obtained via calibration.
Once both the camera and the projector are calibrated, the intersection (or smallest distance) between a vector � ⃗ P j and a vector �� ⃗ C k indicates the local height of the illuminated surface. As a first test, we reconstructed the position of the horizontal plate used for calibration at the different heights. Figure 7 shows the RMS error of the measured values of the calibration plate at different heights when compared to their true values. The average RMS error associated to the technique in this particular configuration is 0.4 mm.

Refraction
In the calibration procedure, the measured surface, camera and projector are in the same medium (air). However, measurements can still be performed when another transparent medium is present without changing the calibration procedure or making a different calibration, for example, when the surface to be measured is underwater. This is  done by accounting for refraction. This change of medium means that the direction of the calculated vectors � ⃗ P and �� ⃗ C will be modified.
To calculate the change of each vector when passing from a medium with a refractive index n 1 toward the medium with the refractive index n 2 , the vectorial form of Snell's law is used, where � ⃗ n is the unit vector normal to the interface between the two media, and � ⃗ v i , � ⃗ v t are the unit directional vectors of the incident and transmitted line in space, respectively. Each incident vector � ⃗ v i , is calculated by defining the direction vectors of both the camera and projector at the interface height H I . This is done using � ⃗ P and �� ⃗ C to locate the position of the crossing (x I , y I , H I ) for each direction vector.
To calculate � ⃗ v t , we express it as a function of � ⃗ v i and � ⃗ n , since these directional vectors lie on the same plane (plane of incidence). This yields where and are unknown. To determine , we take the cross-product of Eq. (2) with � ⃗ n and use Eq. (1), to obtain Additionally, since the directional vectors ( � ⃗ n, � ⃗ v i , � ⃗ v t ) are normalized, the inner product of Eq. (2) with itself yields the following quadratic equation: the solution of which is used to determine . If no solution is found, total internal reflection occurs. By substituting Eq. (3) and the solution of Eq. (4) into Eq. (2), the transmitted directional vector � ⃗ v t for each line in space (for each vector � ⃗ P j and �� ⃗ C k ) is determined. By finding the intersection, or the smallest distance, between the new directional vectors (for both camera and projected dots) the local height of the illuminated surface submerged in the medium with refractive index n 2 is measured.

Reconstruction of objects
To show the capability of the developed technique, the height and shape of different objects and surfaces were measured. The examples treated here show the precision of the method, and for some extreme cases, the limitations.

Experimental setup: rotating table
The experimental setup (shown in Fig. 8) consisted of a 1 × 1 × 0.35 m 3 square Plexiglas tank which was placed on top of a computer-controlled rotating table. The tank was filled to a height H = 21.6 cm with water, with a refractive index n 2 = 1.333 . The position of the interface H I is the height of the water layer, which is here a constant: H I = H . A MEGAPLUS Es 2020 camera, with a resolution of 1600 × 1200 pixels, was positioned at a height H C = 2.06 m above the tank looking vertically down. This gave a field of view of approximately 0.76 × 0.61 m 2 , which was displaced ( 30 cm , 24 cm ) from the geometrical center of the tank. A TH683 BENQ digital projector, with a maximum resolution of 1920 × 1080 pixels, was placed at a height H P ∼ 109 cm from the bottom of the tank and approximately 55.5 cm away from the center of the camera's field of view. The projector was set to an angle of p ∼ 21 • , illuminating an area of 0.73 × 0.375 m 2 inside the field of view of the camera. This defined the effective measurement area or region of interest. Projected images consisted of a series of scattered dots, as shown in Fig. 1. In total, 120 different images were projected for each measurement to densely cover the entire field of view. Two different resolutions were used one with about 50 points per images and one with 300 points per images. These images were generated by creating ellipses with their centroids defined by a uniform grid in the x direction, while in the y direction its position is given by the semi-random MATLAB function "rand". The position of the camera and projector was the same as in Sects. 2.1 and 2.2, giving a spatial resolution of 5 mm and a height RMS error equal to 0.4 mm . This configuration was chosen as default. The effect of the angle of the projector and the size of the projected dots is discussed in Sect. 3.2.1.

Results
Three different types of tests were performed to determine the effectiveness of the measurement technique. First, a well-defined and previously measured surface was chosen to verify that the method was capable of measuring the correct local height. Next, different objects with known dimensions were measured to test the technique under extreme slope conditions. Finally, after the method is shown to work, a sediment bed with different height distributions was measured. For these measurements, the water surface height remained flat at a constant level.
The first surface measured was an undulated plate with a sinusoidal shape. To obtain a reference, the undulated plate characteristics were acquired by taking a calibrated side picture of the plate, obtaining an average crest height of 1.9 cm , a trough height of 1.2 cm and a wavelength ≈ 3 cm . The undulated plate was placed on top of the calibration plate and a false bottom, shifting the surface upward by 1.55 cm . This displacement was measured with a micrometer. This means that the average crest height was h M = 3.45 cm and the trough height h m = 2.75 cm when measured from the bottom of the tank. An image showing the height measurement of the undulated plate is shown in Fig. 9. The measurement yields an average crest height of 3.5 cm , an average trough height of 2.6 cm and a wavelength R = 2.95 cm . It is concluded that there is an excellent agreement both visually and quantitatively with a root mean square error of E RMS ∼ 0.3 mm.

Effects of array configuration and pattern characteristics
Further tests of the measurement technique were done by changing the position of the projector. Three distinct angles were chosen to show how the components of � ⃗ P change. To measure (approximately) the same region of interest, the height of the digital projector had to be changed accordingly. As shown in Fig. 10a, there is an increase in the magnitude in the horizontal components of the vector � ⃗ P = [dx p , dy p , dz p ] for the largest angle . This means that the sensitivity of the measurement technique to changes in height depends on the angle of the projector. Due to the increased sensitivity, a larger displacement in the dot positions can lead to a reduction of the overall error, as shown in Fig. 10b. However, for high , large height variations (or gradients) can cause an overlap of dots when viewed from the overhead camera, significant deformations of the dots, or non-illuminated regions in the surface due to the casting of shadows. This has a negative effect on the measurements, since it can lead to problems identifying the projected dots, an incorrect relation of dots (and positions) for a given pattern, and cause a lack of measurements in the shadowed area. An example of where shadowing and an incorrect dot relation occurs due to missing dots is shown in Fig. 11a where we show an attempt to reconstruct three cylinders of different heights. Although it is possible to mitigate the incorrect identification of dots using a more advanced detection/relation algorithm or by manually relating only visible dots as shown in Fig. 11b, the shadowing cannot be reduced unless the projector angle is decreased. The deformation of the dots can be slightly avoided by projecting smaller dots, but depending on the distortion produced, a different detection algorithm is recommended. Here, the centroids of the light patches in the recordings are found by fitting circles when the distortion is small. Ellipses can also be fitted to the light patches for larger distortions.
Changes in the position and angle of the camera are expected to have similar consequences as changing the position and view of the digital projector. Therefore, it is also possible to increase the sensitivity of the measurements by changing the angle of the camera. However, analogously to the shadowing, there can be a loss of information due to the blocking of the field of view, as shown to the left of the cylinder at x ∼ 150 mm in Fig. 11b. Because of this, the preferred configuration of the camera in this setup is to be looking vertically down and with sufficient height to observe the region of interest. Despite the similarities with the projector, the only way to increase the resolution of the measurement is to use a high-resolution camera or reduce the distance between the camera and the region of interest.
We also examined the effect of the projected pattern, in particular, the dot size on the measurement error. For this, the pixel diameter of the projected dots (as submitted to the projector) was changed from = 4 px to = 40 px while maintaining the distance between them constant ( e = 44 px). The mean diameter of the dots in the images taken by the camera is almost equivalent to those submitted: they range from 6.6 px to 46 px as computed by the detection algorithm. Here, we used the high-resolution configuration with 300 dots per images to ensure that the error is due to the dot size and not to a lack of resolution. As shown in Fig. 12, the RMS error calculated with the mean subtracted remains relatively constant for ≲ 20 pixels but shows an increase for larger values of . The error increases with from ∼ 25 , which is equivalent to a dot size of 30 px in the images obtained by the camera. Since the wavelength of the undulated plate is of about 60 pixels in those images, . 9 a Results of the measurement of an undulated plate, where every dot corresponds to a measurement point from the 120 low resolution images. b Three-dimensional reconstruction of the plate obtained through cubic interpolation. c Comparison of the undulated plate as obtained from a side photograph with the data for 30 < y < 35 mm for low-and high-resolution measurements we conclude that the error increases drastically when the dot diameter is larger than ∕2 . In such cases, the dot size acts as a filter. An additional problem when projecting large dots is that the total amount of light increases, which can lead to detection problems due to, for example, reflections.
The RMS error without subtracting the mean increases for both small and large projected dots. For the small dots, this problem is due to an underestimation of the displacement of each dot's centroid, leading to an underestimation of the measured surface height. Notice that the measurements with = 10 px give the best results. However, the calibration was only performed with = 10 px, probably favoring this dot size. Furthermore, the method to find the centroid by fitting a circle or an ellipse has problems as the dot diameter approaches four pixels, when dots are distorted due to large height gradients, and when there is a large variation in the dots size within a single image due to large changes in height. Fig. 10 a Root mean square of the horizontal components of � ⃗ P = [dx p , dy p , dz p ] as a function of the angle of the digital projector . b Root mean square error E RMS as a function of Fig. 11 Three-dimensional reconstruction of three cylinders obtained when using different methods to relate the dots found in images to their respective lines in space. a Dots in images are related to the closest line in space (after deformation). b Dots are related to their respective line in space by hand selection. The smallest cylinder has a height of 1.1 cm and a diameter of 11.85 cm ; the tallest cylinder has a diameter of 6.25 cm and a height of 18.1 cm while the last cylinder has a height of 17 cm and a diameter of 9.1 cm . Two measurements with different locations for the cylinder are used to show the different types of error Additional changes can be made to each pattern by increasing the separation e between the dots. This effectively changes the resolution per image. The resolution can be increased as long as dots do not overlap or their position is not inverted. For these cases, a larger distance between dots (a smaller density of dots per image) can improve the measurement since the possible locations to which a dot can be related are reduced. Reducing the number of points per images can be compensated by increasing the number of images.
While the values of the error as a function of the dot size and the projector angle are limited to the setup described in Sect. 3.1, it is possible to apply the same approach regardless of the configuration. Therefore, it is recommended that the angle of the projector should be set based on an initial estimation of the height distribution of the measured surface. This way, the sensitivity of the method to changes in height is optimized for the setup while avoiding the casting of shadows. The projected pattern characteristics should be chosen based on the camera position (i.e., the resolution) and the dot size relative to the typical horizontal length scales to be measured. Depending on the pattern chosen and its distribution, a number of projected images can be set to determine the horizontal resolution of the measurement. The lower the density of the dots per image, the larger the number of images that is needed to maintain the horizontal resolution. However, there might be restrictions in the amount of images that can be projected in a given time due to both the acquisition and projection rates and typical time scales of the flow above the measured surface.

Granular surface measurements
To test the measurement technique on a granular surface, small white polystyrene particles were added at the bottom of the experimental setup until a sediment bed of thickness h ≃ 1.6 cm was formed. The particles have an average diameter D = 280 μm and a density p = 1.2 g cm −3 . Due to the density difference, particles settle and produce a sharp interface with a well-defined surface.
Since large-and small-scale structures can develop in a sediment bed, two different particle bed distributions were tested. To determine if small-scale structures can be properly measured, a "flat" particle bed was measured before and after being disturbed by a flow over it. The reconstruction of the sediment bed is shown in Fig. 13. As can be seen, the technique is capable of measuring the ripples produced after the disturbance, which had a height of less than 1 mm.

Reconstruction with a curved water surface
To demonstrate the possibility of measuring through an interface of any shape, a test surface was measured under a parabolic interface. To do this, the experimental setup was set to rotate at a constant angular velocity . As described in the classical problem known as "Newton's bucket", the free water surface of a rotating container follows a parabolic curve given by: where g is the gravitational acceleration, r = √ x 2 + y 2 is the distance from the axis of rotation, and H o is the height of the water layer at the center when in rotation. This height is determined by volume conservation since the same volume of water is contained in the tank with or without rotation. For a square tank with side L, this means that L 2 , Fig. 13 Three-dimensional reconstruction of an initially "flat" sediment bed a before and b after a producing a local disturbance. Both images show the reconstruction of features even with small heights in the order of 1 mm . The increased height measured at x = 0 is due to the presence of a plate used for reference where H is the height of the water with the tank at rest. Using Eqs. (5) and (6), the water surface height is calculated throughout the domain. The water surface height in the measurement region is shown in Fig. 14. Implementing this position-dependent height in the refraction calculation of each line in space allows measuring submerged surfaces inside the rotating tank, for any value of . To test this approach, the undulated plate is used. The initial water layer height was the same as in previous measurements, H = 21.6 cm . To produce a pronounced parabolic deformation of the free water surface, as set to ∼ 1.5 rad s −1 , leading to H o = 19.7 cm. Figure 15 shows the reconstruction of a section of the undulated plate under the parabolic water surface, together with the reconstruction of the plate under a flat water surface. Both reconstructions have a shape similar to the actual surface, but they show an almost constant height difference. The height of the undulated plate measured under a parabolic water surface is underestimated throughout the domain by 2.5 mm . While the underestimation diminished when reducing , the correct height was only recovered when = 0 rad s −1 . When removing the average height before computing the root mean square error E RMS of the reconstructed surface under rotation, the error decreased considerably to ∼ 0.5 mm , as shown in Fig. 16. This indicates that the underestimation in height is systematic.
By artificially changing the value of H o obtained with Eq. 6, we found that the RMS difference (without removing the averaged height) decreases to less than 1 mm when H o ≈ 19.25 cm , as shown in Fig. 17. As H o increases or decreases from this value, the overall error of the reconstruction increases linearly with a slope of about 0.475. Meaning that an error of 1 mm in the position of the interface results in an error of less than half of that in the measurement. This effect of the error in the location of the interface on the accuracy of the measurement depends on the total depth (as also found by Akutina et al. 2018, for 3D particle tracking velocimetry through an interface) and the angle of the projector. A uniform error in the water height results in a tilt in the reconstruction. However, this tilt can be easily corrected for if there is a reference height.
The error in the results obtained using the value of H o calculated with Eq. (6) is probably due to an error in the computation of the total water volume and the rotation rate of the table. In particular, we observed that the side walls of

Measurements through a time-varying water surface
It was shown in the previous sections that the photogrammetric technique is capable of measuring the height of submerged surfaces under a static water layer as long as the interface height is known. However, it is also of interest to consider how to implement this technique in cases where the submerged surface to be measured, the interface, or both vary in time.
An example of interest is the evolution of a sediment bed under an oscillating channel flow, where both the interface and sediment bed vary in time. Here, a sinusoidally varying water height with a period T forces a flow through a channel. This flow exerts a time-dependent shear stress (t) on the sediment bed. If the shear stress exceeds a critical value (t) > cr , particles are set into motion, leading to the local erosion and accumulation of sediment across the channel.
To test the developed measurement technique with an oscillating water height and an evolving sediment bed, a second experimental setup was used. In this setup, the photogrammetric technique was combined with PIV to measure both the flow velocity field and the evolving sediment bed height distribution. Figure 18a shows the setup used to test the photogrammetric technique under a moving water surface. It consisted of a 1.5 × 1.0 × 0.3 m 3 perspex tank, partially filled with water to a height H. Inside the tank, a plate of 1.0 × 1.0 × 0.08 m 3 was placed as a false bottom. The water layer thickness was between 5 and 8.5 cm above this false bottom. Two curved metal barriers were positioned on top of the false bottom to form the tidal channel with a length of 65 cm and width of 21 cm . This defines our region of interest. The channel openings were rounded, with a radius of 9 cm , to avoid erosion on the sediment bed due to strong convergence and separation of the flow. The tank was divided by the channel into an open area, representing the sea side, and a semi-enclosed basin or estuary.

Experimental setup: tidal channel
Inside the channel, an oscillating flow was produced by inducing a change in the water level at the sea side with a vertically moving piston. The piston was partially submerged 3 cm in the water and moved sinusoidally up and down with a period T = 21 s and amplitude L piston = 3 cm . The cross-sectional area of the piston A piston = 1000 cm 2 remained constant throughout all experiments. Two flow reversals occur during each period, transitioning from ebb to flood at t = (2n − 1)T∕2 and flood to ebb at t = nT with n = 1, 2, 3, ... Over the false bottom, opaque cylindrical polystyrene particles were added until a 2 cm sediment layer was formed. Particles have a density of p = 1.055 g cm −3 , with particle size distributions D 10 = 1.7 mm , D 50 = 2.1 mm , D 90 = 2.9 mm , and maximum size of D 10 = 4.0 mm . The small density difference between water and sediment was chosen to increase the mobility of the particles, making the sediment transport easily noticeable.
On top of the setup, two cameras were placed at a height H C = 177 cm above the false bottom, looking vertically down. Camera 1 (ALLIED PIKE) has a resolution of 1388 × 1038 pixels, while camera 2 (MEGAPLUS Es 2020) has a resolution of 1600 × 1200 pixels. Both operate at 30 fps. A MH530 BENQ digital projector operating at 60 fps was mounted behind the tank at a height H P = 161 cm . Since the positioning of both camera and projector was different from the experimental setup described in Sect. 3.1, new calibrations were made for both camera and projector.
Measurements of the flow velocity were performed using PIV. For this, the water surface was used as the plane of measurement to avoid interference from suspended sediment particles. As tracer particles, we used COSPHERIC UVY-GPMS-0.97 600 ∼ 710 μm low density, fluorescent, polyethylene micro-spheres, with a density t = 0.980 g cm −3 . These particles are illuminated by four pulsed UV LED lights (Lz4-44UV00, Led Engin) which are positioned next to the channel (as shown in Fig. 18b). Ultraviolet illumination was used to cause the particles to emit a yellowgreen light (with wavelength 365 nm ). Calculations of the velocity fields are performed using commercial software (PIVview3C), developed by PIVTEC.

Morphodynamic measurement implementation
To implement the photogrammetric technique on this timevarying problem, the projection and recording time must be small in comparison with time scales of the flow and the sediment bed evolution. This can be achieved by reducing the number of projected dot patterns or increasing the projection rate. In this setup, to reduce the time required for measurements, only 15 dot patterns were projected. These patterns consisted of an array of 10 × 20 dots arranged in a square grid, with a dot size (at the sediment surface) of 2.5D 50 . Each consecutive image is slightly shifted in the x-and y-direction to measure over the entire channel. This means that the region of interest is covered by 3000 measurement points, giving a spatial horizontal resolution of 4.5 mm . Figure 19 shows a photograph of the first projected image with a schematic of the change in position for each consecutive image. The order of projections was chosen to average out a small systematic error in the measurements produced by the changing water level observed when testing the setup. Images are projected/recorded during each flow reversal, where v flow ∼ 0 cm (causing (t) < cr ) to avoid movement in the particle bed while performing measurements. Due to a slow communication between the program submitting the images and the projector, the working frame rate of the projector (60 fps) could not be used for the measurements. Each dot pattern was projected and recorded within 0.075 s meaning that the total measurement time was 2.1 s during each flow reversal (i.e., 10% of the oscillation period). However, this time could be reduced using better hardware and a more efficient way of submitting the images. The change in the water level height H , is given by the ratio between the water volume displaced V = A piston ⋅ L piston and the total area of the domain A tot , such that: The change in height propagates along the channel with a wave velocity c w = √ g ⋅ H W ≈ 100 cm s −1 , as obtained using a shallow-water approximation. This means that the time required for a surface disturbance to travel across the full domain is approximately 2 s . Since this time is much smaller than the oscillation period, it is assumed that the water level rises and decreases simultaneously across the setup. This means that measurements can be taken under an approximately flat layer of water. However, due to the flow, small disturbances in the water surface are still expected.
To determine the effect of the changing water level on the measurements, the height of the bottom ( h 0 = 0 cm ) is computed with a stationary water level H = 6.3 cm and during flow reversal with two different mean water levels H = 6.3 and 7.5 cm. As seen in Fig. 20, the measurements taken during flow reversal differ from those under the stationary water layer. During flow reversal, an offset of 0.5 mm is observed over the entire length of the channel, with an increase in the scatter of the calculated height. This increase in scatter can be explained by the small disturbances at the water surface, while the small offset can be explained by the small increase in the overall water level. A similar effect is observed in Sect. 4, where there is an underestimation of the overall height measured caused by an incorrect water level height. Since the offset and slight increase in scatter are small compared to the initial bed configuration (differing by almost two orders of magnitude), the error is considered within the bounds of the measurement technique.

Implementation of PIV measurements
It is now common to use PIV to measure a two-dimensional velocity field of a flow. Hence, only a brief description including the most relevant aspects is included here; for more detailed information, see Raffel et al. (2007). PIV consists of capturing a moving distribution of tracer particles as a sequence of images. Each of these images is then subdivided into smaller windows which are correlated with the ones in subsequent images taken a time period t later. Here, the minimum value of t is determined as the frame rate of the camera. However, additional restrictions are imposed because measurements of the sediment bed also take place in the same region of interest. Ideally both measurements should take place simultaneously. For this, the illumination and capture of images for one technique should not interfere with those of the other. To remedy this, images are captured out of phase with those required by the other measurement technique. In essence, this means that the projected dot patterns used for the morphology measurements are recorded only when the tracer particles are not illuminated and vice versa.
To acquire the morphology and PIV measurements out of phase within a very short time span, we use the fact that the digital projector does not illuminate continuously. In fact, it projects twice the same image every 1∕60 s with a time interval t p between each projection. In this way, recordings for the morphology measurements are taken when the patterns are projected, while images for the PIV technique are taken during the time interval t p when no light is emitted by the projector. To synchronize the projector, UV LED lights and cameras, a function generator and a triggering box are used. Both cameras and the power source of the UV LED lights are connected to the trigger box, while the trigger box itself and the digital projector are connected to the function generator. For these experiments, the digital projector is used to synchronize both cameras (operated in asynchronous edgetriggered mode) and the LEDs by sending a signal (pulse) when the projection image changes. This signal is given by the vertical sync (Vsync) of the projector and is related to its frame rate. This type of interlinked image acquisition has also been applied by Dalziel et al. (2007) for simultaneous PIV and synthetic schlieren measurements, where an LCD monitor was used instead of a digital projector to synchronize the cameras. By choosing to project only secondary colors (like yellow), the time interval t p is extended slightly from 3 to 4.5 ms. The activation by the Vsync occurs as follows. When a Vsync is detected, the function generator sends two block pulses to the trigger box, one directed to the camera used for PIV (camera 2) and the second to the power supply of the UV LED lights. The UV LED lights are activated for 4.7 ms , overlapping the exposure time of the PIV camera ( 4.5 ms ). The camera used for morphodynamic measurements is triggered after the UV LED lights are deactivated, with an exposure time of 15 ms to capture the two consecutive projections of the same dot pattern. The pulse for the cameras is extended beyond the next Vsync to force the system to a constant 30 fps. The activation of the cameras occurs only when the trigger signal changes from being under to being over 3 V. Both the PIV camera and the UV LED lights are active when the digital projector is not projecting and image. The camera used for morphodynamic measurements captures two consecutive projections of the same dot pattern.

Results
For the experiment, a flat 2 cm sediment layer was made in the channel. Disturbances are then produced in the bed to observed their evolution in time. Measurements of the sediment bed in the channel are taken for over 500 tidal periods. Two measurements of the sediment bed are performed after every five periods during the transition from flood to ebb and during the transition from ebb to flood. An example of an initial disturbed sediment bed height h(x, y, t = 0 s) with H = 7.6 cm is shown in Fig. 21. In this figure, each dot represents a height measurement, showing the location of three heaps of sediment: a small one near the entrance of the channel to the semi-enclosed basin ( x ∼ 85 cm ) with a height h ∼ 3.5 cm and two large ones at x ∼ 110 cm and x ∼ 130 cm with maximum heights of 4.2 cm. Figure 22 shows the evolution of the sediment bed. For the case shown, the patch initially located at x ∼ 110 cm migrated in the negative x-direction. Furthermore, it grew as it traveled, filling the complete width of the channel within 150 tidal periods. The migration speed of this patch was approximately 6.5 × 10 −5 cm s −1 or 0.014 cm∕T . A similar behavior was observed for the patch initially at x ∼ 130 cm , growing in height and width as it moved toward the left.
The vorticity fields obtained using PIV for t = 5 T, 125 T, 250 T , and 375 T are shown in Fig. 23. An increase in vorticity can be seen in a region near x = 113 cm . In the beginning of the experiment t = 0 T, 5 T , this area corresponded to the maximum height in the sediment bed. Once t > 250 T , the highest vorticity remained near x = 110 cm where the height of the sediment bed is the lowest. Measurements of the water surface velocity field were performed at t = 0 T, 5 T, 125 T, 250 T, 375 T and 500 T . The velocity averaged across the channel at x = 105 cm is shown in Fig. 24, where an asymmetry between the magnitude and duration of the flood and ebb phase is observed. The magnitude of the velocity during flood is 2.8 cm s −1 , while during ebb, it is − 3.0 cm s −1 . The duration of the flood and ebb was 11 and 10 s, respectively.

Discussion
As discussed by Porter et al. (2014), the main drawbacks of using photogrammetry to measure a sediment bed under a water layer is the temporal resolution and the intrusiveness of the cameras on the flow or the refraction caused by the air-water interface. However, with the improvements made to the photogrammetric technique, these issues can be overcome. This can make the photogrammmetric technique presented here a viable measurement method for submerged, evolving surfaces, as long as the measurements occur in shorter time scales than the changes in the surface to be measured.
As shown in previous sections, the developed photogrammetric technique is capable of reconstructing a variety of surfaces even while submerged in water. Only one calibration for the camera and projector is needed to perform experiments with different water levels as long as the camera and projector are kept in the same position. For the setups used in this study, the RMS error was less than 1 mm . This demonstrates the robustness of the technique in dealing with refraction at an interface. This is relevant for experimental morphology measurements, since cameras do not need to be placed underwater disturbing the flow or the water does not need to be drained to perform the measurement.
However, if the location of the interface is incorrect, the calculated components of the refracted vectors will be affected, leading to an incorrect height reconstruction. For the photogrammetric technique, the error depends on both the water depth and the projection angle. In the current paper, we estimated indirectly the position of both a parabolic and a moving water surface. For more accuracy or Fig. 21 Example of the height distribution h(x, y, t = 0 s) of the sediment bed at the beginning of an experiment. The graph shows the 3000 measurement points taken of the surface of the sediment bed more complex varying water surfaces, the technique could be combined with measurements of the water surface using, for example, synthetic schlieren (Moisy et al. 2009).
The choice of camera and projector and their location (including the projection angle), the size and distance between the dots, the location of the interface, the material composing the surface to be measured, the vertical and horizontal size of the disturbances to be measured, they all affect the accuracy of the method. This means that there is a large amount of possible configurations, and hence, a welldefined set of design rules would require further research. However, some guidelines can be derived from the present paper: (1) a larger angle of the projector with respect to the vertical results in a more sensitive measurement, but it finds limitations due to the reflection of light from the interface and the possible casting of shadows from the patterns that emerge, and (2) the diameter of the dots should be at least one-third of the horizontal length scales to be measured, about three times larger than the typical grain size so as to filter those length scales, and larger than about five pixels (both in the images sent to the projector and those captured by the camera) to avoid detection problems.
In a way, the fact that we cannot provide a complete set of design rules results from the flexibility offered by the method, particularly, due to the use of a digital projector. This flexibility is extremely helpful since the size of the dots, their separation and their distribution can be optimized for each surface measured. In addition, multiple and different pattern configurations can be used without any difficulty. For example, in the current paper, we used two types of patterns: random and ordered ellipses. Ordered patterns can achieve a larger dot density and their characteristics are easier to define, while a random distribution can improve correct detection and reduce aliasing when the apparent movement of the dots exceeds the spacing given by the projected pattern (as also discussed by Dalziel et al. 2007, for a pattern Fig. 22 Bed morphology measured during the transition from flood to ebb at t = 5 T, 125 T, 250 T and 375 T . It is observed that the initial perturbations migrate and grow, covering the entire width of the channel after 150 periods matching technique). However, there does not seem to be, in general, a clear superiority of either one. What is true for both is that the configuration of the pattern must be made with care. This is because, depending on the height distribution of the measured surface, projected dots can merge or vanish due to shadowing. If not accounted for, the lack of dots (or overlapping) can lead to non-measurable areas or an incorrect correlation of a detected dot to its originating ray.
Furthermore, using a digital projector allows for the projection of many patterns in a short time. For example, the projectors mentioned in the current paper could conceivably run at 60 fps. This allows for the measurement of a nonstatic sediment bed, provided movement occurs due to bedload transport, that there are clear-water conditions, that the cameras work at the same frame rate, and that the changes in the sediment bed and the free surface are relatively slow compared to that frame rate. The projection of many patterns could also allow to project patterns with multiple characteristics (e.g., dot size) and chose, a posteriori, the best one or even a combination of patterns. Digital projectors also allow for the projection of colored patterns, which were not discussed here, but Kolaas et al. (2018) showed that monochromatic light can improve the measurements of a water surface since light diffraction is wavelength dependent. In our first experimental configuration, no significant difference was observed in terms accuracy in preliminary experiments. It is thought, however, that this was due to the type of camera or to the colors chosen to maximize the brightness of the dots, which reduced the displacement. However, it is expected that the light wavelength will have a larger effect for monochromatic patterns (especially those in the extremes of the visible light spectra) or for larger water depths.
The use of a digital projector comes with a caveat. Since projectors are not designed for precision measurements, certain features that might improve their intended function as projectors for presentations and movies might interfere with the accuracy of the measurements. In particular, we observed that the projected images were, sometimes, changing in time. It is then important to first test the behavior of the digital projector.

Conclusion
This paper describes the implementation of the photogrammetric technique using a single camera and a digital projector. Improvements on the original method are performed by reconstructing the light rays that emerge from the camera and the projector and using the intersection of these lines in space to measure the local height of a surface. The developed method also accounts for refraction effects when a different medium is present using Snell's law to calculate the deviation of the light rays at the interface between media. Deformations due to the projection angle and due to the camera lens are also included by performing calibrations of the digital projector and camera.
Since a local implementation of Snell's law is used, a correct reconstruction of the light rays from both projector and camera is obtained even through non-flat and varying water surface. This is shown by performing measurements through a parabolic and a vertically oscillating air-water interface. Additionally, both the morphology and velocity fields of the flow can be obtained by alternating/sequentially capturing images for the photogrammetric technique and for PIV.
The error of the measurement technique for the configurations used in the current paper was E RMS < 1 mm even when working with sediment particles. However, the accuracy should be estimated for each setup due to the large flexibility offered by the technique. In fact, the accuracy depends on several factors, such as the angle of the projector, the type of surface measured, the camera resolution and location, the chosen pattern, the water height, and the accuracy in the location of the interface. Although the dependence of the error in the measurement as a function of some of these factors was explored, the particular values are limited to this setup. While more research is needed to give general rules for an optimum implementation, this paper illustrates what the technique can achieve.
As shown, the developed photogrammetic technique produces accurate, non-intrusive, high spatial and temporal resolution measurements of a surface, even when submerged in water. Furthermore, the method allows for a tailored approach for its implementation in an experimental setup, making this a valuable technique for the measurement of evolving sediment beds.