Concepts for compensation of wave effects when measuring through water surfaces in photogrammetric applications

A common problem when imaging and measuring through moving water surfaces is the quasi-random refraction caused by waves. The article presents two strategies to overcome this problem by lowering the complexity down to a planer air/water interface problem. In general, the methods assume that the shape of the water surface changes randomly over time and that the water surface moves around an idle-state (calm planar water surface). Thus, moments at which the surface normal is orientated vertically should occur more frequently than others should. By analysing a sequence of images taken from a stable camera position these moments could be identified – this can be done in the image or object space. It will be shown, that a simple median filtering of grey values in each pixel position can provide a corrected image freed from wave and glint effects. This should have the geometry of an image taken through calm water surface. However, in case of multi camera setups, the problem can be analysed in object space. By tracking homological underwater features, sets of image rays hitting accidently horizontal orientated water surface areas can be identified. Both methods are described in depth and evaluated on real and simulated data.


Introduction
The article addresses a common problem when imaging and measuring through moving water surfaces.Due to the refractive properties of water, standard photogrammetry methods cannot be applied without modifications.While simple approaches try to compensate effects of refraction by introducing empirical parameters (Woodget et al., 2015;Partama et al., 2017), more sophisticated methods work with strict geometric modelling of refraction (Mulsow & Maas, 2014;Maas, 2015).The latter ones require information about the shape of the water surface at the time of image capture.In the simplest case, the water surface can be described as a horizontal plane.This is sufficient if the water surface is calm (Mulsow, 2018).Nevertheless, many measurement tasks, both in nature and in laboratory environments, have to be carried out with a wavy water surface.Furthermore, the shape of the water surface is instationary and constantly changing over time.In the best case, the water surface can be modelled as a parametrized surface with a limited number of parameters (e.g.horizontal plane with time-dependent change of water level, moving sinusoidal wave).The problem may become much more challenging with real world water surfaces of rivers, lakes and the sea.Complex water surfaces can usually not be described mathematically by a global parametrization.Therefore, methods based on a strict global modelling (e.g.Fourier surfaces) of the water surface cannot be applied in such cases.A strategy to overcome this problem is to switch to a local modelling of the water surface (Yiming et. al. 2018).This relies on local smoothness and local low complexity.An Achilles Heel of such methods is the low redundancy in parameter determination -only the image observations through a small portion of the water surface can be used for the reconstruction of the water surface, which introduces a number of additional unknown parameters in addition to camera orientation parameters and water bottom point coordinates.
In this article, we present strategies for the compensation of wave effects when imaging objects through a non-planar nonstationary water surface.Generally, our methods base on short image sequences taken from a stable position rather than on single images.Obviously, this can only be realized in certain applications (e.g. with static camera setups or with multiple hovering (i.e.stationary) drones), but then the processing of short image sequences offers a way to avoid the challenging modelling of the water surface.The presented compensation methods can be categorized in two main sections: Data analysis in image (Section 2) or object space (Section 3).

Analysis in Image Space
As shown in Figure 1, a non-planar non-stationary water surface causes disturbing effects when imaging submerged structures.The most obvious one, at least in this example, is the glint effect from specular reflection of light.Another one is the water bottom distortion effect caused by refraction while imaging through turbid water -water bottom points will move over time in the image of a stationary camera.Such freeform water surfaces are difficult to model, making image correction based on a single image and a global surface model a rather complex task.By analysing a short image sequence taken from a stable position, a major improvement of the imaging quality can be achieved.After the correction procedure described in the following, a glint-free image with a quality and geometry (almost) as if acquired through a planar water surface can be provided.This can then be forwarded to standard photogrammetric procedures, such as twomedia forward intersection.

Method
Our method is based on the assumption that the shape of the water surface changes randomly over time.Therefore, for a specific pixel position at least one image of s sufficiently long image sequence should provide the "correct" (glint-free, distortion-free) grey value.Further, we use the fact, that the water surface moves around an idle-state (calm planar water surface).Thus, moments at which the surface normal is orientated vertically should occur more frequently than others.The simplest method to correct for the disturbing effects is to apply median-filtering on the time-axis of the grey-values in each pixel position for all images of the sequence.Values of 255 (overexposure) are excluded from median-filtering.For pixels showing mainly overexposed values, the algorithm has to be adapted.In that case, the Minstore-value (smallest pixel value) can be taken for the corrected image for example.If only fully saturated values of 255 were captured for one pixel position, the corrected value is interpolated from neighbouring pixels.First tests have proven the concept (see Section 2.1), and the method was transferred to a practical application (see Section 2.2).The image content corresponds to an image of a projection through a calm planar water surface.

Evaluation
For the evaluation of the method, a dataset recorded by the Department of Geodesy and Geoinformation (Vienna, Austria) was used.As shown in Figure 3, the water body is simulated by a tub (mortar pan) filled with two layers of stone and water.Two different camera setups were realised -configuration A: 4 synchronized grayscale cameras in a 2x2 nadir arrangement; configuration B: 2 nadir and 2 oblique cameras.The cameras were triggered simultaneously by an Arduino mini-computer every one second.Water waves can be induced repeatably by a small water pump located on the edge of the tub.First, an image set was taken simultaneously through a calm water surface as a reference.Further, a second data set was acquired through wavy water from the same camera positions.
Figure 7 shows the range of displacement of specific features inside an image sequence of a nadir camera.So, when using the image coordinates of a feature from a single (wavy) image for simple forward intersection assuming a flat water surface, a significant error for the position estimation in object space can be assumed.In this specific case, for a stereo pair with an imaging distance of 1.6m, a base of 0.4m, a focal length of 8mm and a pixel size of 3.45µm, an accuracy in height of 1.4mm can be estimated for a horizontal still (planar)water surface (assumed accuracy for image measurements: 0.5pixel).However, an imaging error of 5 pixel in the stereo images would mean a loss of height accuracy of factor 10 (14mm).
The images from one nadir camera (moving water) were processed according to the procedure as described in the previous section.For evaluation, the derived image (corrected) was compared with the image taken through calm water (reference).
As comparison criteria, the relative displacement of image features was used.Interest points were extracted in the reference image first.Then these points were matched in the corrected image via Least Squares Matching (LSM).An RMS of 0.4/0.5pixel in x/y position could be achieved, with the residuals showing a normal distribution (see Figure 9).When comparing the RMS with the range of movement inside the image sequence (standard deviation of ~5pixel, see figure 7), one can see that the effect of waves could be reduced by a factor 8-10.This means an increase of accuracy for the height component of ground points by the same factor.
In addition to the increase of accuracy, further advantages of the method can be identified.The correction can improve the imaging quality and completeness compared to single images (see Figures 1 and 2), as glint effects can be effectively eliminated.Further, the use of images freed from wave effects can ease the orientation and 3D model extraction.The whole modelling of refraction can be reduced to a planar air/water interface problem.Therefore, when only the geometry of objects (in our case the river bed e.g.rocks) has to be determined, the complex modelling of the spatio-temporal wave movements can be omitted.Additionally, the image sequences don't have to be acquired simultaneously.Therefore, image data can also be acquired by taking short images sequentially at multiple positions, thus reducing the complexity (no triggering) and costs for equipment.

Practical Application
Multi-camera systems have been successfully used in hydraulic engineering experiments to measure riverbeds (Radice &, Zanchi, 2018).For instance, the Federal Institute for Hydraulic Engineering (BAW) in Karlsruhe has been using an optical system (AICON/Hexagon) to measure the riverbeds of hydraulic models for many years (Godding et al., 2003).The measuring principle is essentially based on a pattern projection with simultaneous multi-image recording (3 cameras, See Figure 10).
The system can also measure water bed topography through the open water surface.However, the modelling of refraction in the software of the turnkey photogrammetric measurement system is limited to calm water surfaces (planar horizontal interface).In the case of moving water surfaces, effects such as varying surface normal together with overexposure (glint), smeared pattern projection and motion blurring often lead to a massive loss in accuracy or incomplete waterbody bottom models.Therefore, the water has to be perfectly calm when measuring.This is obviously hampering the simulation of dynamic processes.For each measurement epoch the whole experiment will have to be paused.After stopping the water flow, the water movement has to calm down.This takes up to 30min before the actual measurement can be started.Besides the time issue, the stop-and-go procedure may falsify the whole experiment.
Consequently, there was a need on the user side for a more practical solution.The Institute of Photogrammetry and Remote Sensing at TU Dresden was commissioned by the BAW to find a strategy to adapt the existing system to moving water surfaces.At first, the data acquisition procedure was modified inside the control software from single image to short image sequence recording.The images are forwarded to the image correction procedure as described in Section 2. Then the corrected images (see Figure 2) were fed back to the system software for further processing.Inside the software, the corrected images were treated as if they were standard images taken through calm water.First results proved the concept.The degree of model coverage could be increased significantly.It is noted that a slight loss in accuracy with respect to strict time-resolved geometric water surface modelling still has to be expected.A more detailed evaluation is still work in progress.
Figure 10.Left: 3D water bottom measurement system at BAW in Karlsruhe, Germany -grid projector (P) and three synchronised cameras (C1-3).Right: Projected grid together with a control point.Note the low texture of the gravel material showing the need for structured lighting.

Analysis in Object Space
While the method presented in Section 2 analyses the change in the pixel value over an image sequence, the second (alternative) approach presented in this paper transfers the problem and solution into object space.The basic idea of the approach is to find combinations of images from multi-view image sequences, wherein the water level is horizontal.We assume that at least two short image sequences from different cameras or positions are available.With known interior as well as exterior orientation for each image, corresponding image ray vectors for a chosen ground feature can be projected into object space.These vectors can be traced through a pre-defined plane (horizontal surface of calm water at a well-defined level) towards the water bottom.In a following step, the distances (deviations from ideal intersection) between corresponding vectors (same point in different cameras) are analysed (see also Figure 15).Vector-sets with small distances (which means a small standard deviation for the intersection point) will be considered as candidates for further analysis.This basic idea is applied in a way that tentative intersections between combinations of images of short image sequences of at least two cameras are performed.From these combinations, the best fitting vector set for the centre of each ground feature is handed over for a final spatial intersection for the determination of the 3D-coordinates of that ground feature.The whole process relies on the fact, that the probability of a combination of best fitting image rays is highest for local horizontally orientated surface patches.As well as for the image based method, the analysis in the object space does not require simultaneous captured image sequences.Therefore, only one camera is required for data acquisition.

Evaluation based on simulated data
So far, we analysed the approach as described in this section only on simulated data.The numerical simulation of water waves is a common task in computer graphics.Some widely used methods base on Fast Fourier Transformation (FFT).(Tessendorf, 2001) presented a mathematical toolbox for realistic water body simulation following this principle, which has also been used for analyses of wave effects in lidar bathymetry (Westfeld et al., 2017).In this, the wave height at a specific XY-position is represented by a sum of sinusoids with complex, time-dependent amplitudes.The wave spectra can be derived from theoretical as well as measured data, such as the Phillips spectrum (Tessendorf, 2001).In this, wind speed, wind direction as well as gravitation are taken to account.In order to obtain a realistic appearance of the wave surface, a random 2D-field (Gaussian random numbers with mean of 0 and standard deviation of 1) is introduced in to the FFT.FFT-based random ocean surface models are easy to implement, realistic and real-time capable.However, these models are limited to open ocean scenarios.The interaction between water and obstacles cannot be handled directly (Jeschke & Wojtan, 2017).Further, the height field cannot be parametrized like a surface function (Z = f(X,Y) ) due to the random 2D-field inside the model.For the tests the Ocean Wave Simulator after Tessendorf was implemented.Based on the calculated height field, images of the water surface can be simulated (see Figure 11).However, the calculation of image coordinates of submerged objects is more complex due to refraction.The imaging path from the object space to the image space can be calculated via ray tracing algorithms as presented in (Mulsow & Maas, 2014).However, these algorithms are limited to parameterized surfaces.Due to the non-functional character of the Tessendorf simulation, the ray tracing has to be modified to meet random 2.5D surfaces.Similar to hierarchical matching procedures, the ray tracing was performed sequentially through models of the water surface from low to high resolution.In each resolution step, an 8-parameter Fourier surface was fitted on to the surface points and the ray tracing was performed.The calculated interface point was then used as the starting point for the next step.In the highest resolution the water surface is approximated by a cubic spline surface (9x9 neighbourhood).The whole ray tracing process cause a significant expense in computational load.
For the tests a moving water surface was simulated on a grid of 20x20m with a mean level of 0 m and a min/max variation in height of -0.5/0.5m.Virtual cameras were placed 10m above the mean water level in a quadratic pattern (2.5m x 2.5m distance).
The sensors were introduced as standard machine vision cameras with a focal length of 30mm, a resolution of 4288x2848 pixel and a pixel size of 5.5µm.Figure 14 shows the simulated water surface together with the 4 virtual cameras.
To quantitatively evaluate the method presented in Section 3 (identification of image rays intersecting locally horizontal surface patches), the same setup as mentioned before was used, but now with four cameras (see Figure 12) and just a single underwater point.Water surface states as well as images were simulated over a time span of 100sec in 1sec intervals.Figure 12 shows the dynamics of the water level change for the intersection point of an image ray with the water surface.The water level varies from -0.4 to 0.5m.This results in a variation of the underwater point position in the image up to ±1mm (±180pixel, see Figure 13).First, the number of candidates for finding the best fitting group of image points (i.e.image rays) was reduced by analysing the distance from the mean of image coordinates of the underwater point for each camera.Here, the standard deviation can be used for separation for example.For this experiment the number of candidates could be reduced by 70-80%.The ray tracing is performed for the remaining points from the image space through a plane water surface (horizontal orientated at level of 0m) for each image point coordinate.
In the following step, the distances between image rays of two cameras are calculated and the pair with the smallest distance is taken for forward intersection (first approximation of underwater point coordinate).The distance from that point to the image rays from the remaining cameras are calculated and analysed.The bundle of rays (one for each camera) with the smallest distance to each other are taken for a final integrated forward intersection.Finally, the resulting point coordinate together with the standard deviation can then be compared with the reference.The whole virtual experiment was carried out for various point positions, with error free image coordinates as well as with superimposed normally distributed errors.
In order to give an impression of the accuracy potential of the method, the simulated image measurements through planar water surface were distorted with a normal distributed error of 0.5 pixel.These data were handed over to multi-ray forward intersection and the calculated underwater point coordinate was compared with the reference.Here, an RMS of 1/1/7mm could be reached in X/Y/Z which marks the lower limit of accuracy.Then error free image measurements of the wavy water surface were processed as described before (finding best fitting quadruple of image rays).An average RMS of 4/4/20mm could be reached for the described setup.This means a drop of accuracy by a factor of 4. In order to give a more realistic accuracy estimation, the image coordinates were distorted again with a normal distributed error of 0.5 pixel.This resulted in a further lowering of the accuracy down to 8/8/41 mm.As a rule of thumb, a drop of accuracy of factor 10 against the highest possible accuracy can be given.When comparing the points calculated via multi-ray forward intersection based on uncorrected image measurements (one quadruple for each time step) against the reference, a RMS of 114/144/835 mm could be reached.These drastic values show the need for a proper correction routine, at least in this case.
However, these values are valid for this specific setup with relatively large waves.In the case of more calm waves, the difference should be less but still significant.

Optional refraction-free approach
Another variation of the discussed method would be the search for moments of collinearity between the water surface normal and the incoming image camera ray.In that case, there is no effect of refraction, and the image ray can be seen as straight from the image point to the corresponding object point as for an 'in air' scenario (see Figure 16).However, as the occurrence of such moments in multiple views is much less likely than those of horizontally orientated surface patches, the computational effort of this approach would be rather large.So far, we have not worked on an implementation.
Figure 16 -Principle of the refraction-free approach -the search for moments of collinearity between the water surface normal and the incoming image ray.In that case the refraction has no effect.

Conclusion
In this article we present two different methods for compensating for wave effects when measuring through moving water surfaces.Random wave movements cause refraction effects that are difficult to describe mathematically in space and time.Therefor we concentrated that do not require modelling of the water surface.
The whole problem is reduced to a simple planar water surface as discussed in the previous chapters.Basically, we are searching for entities in the image data material that correspond to a measurement through a flat water surface.The search can be performed in either image (chapter 2) or object space (chapter 3).In simulations, both strategies showed their potential to overcome the massive distortion issues when measuring through moving water surfaces.However, so far only the correction method based on analysis in image space has been evaluated in real world experiments.The big advantage of both methods is that the photogrammetric measurements through moving water can be carried out sequentially.This means that a single camera is sufficient for such tasks.The expense for synchronisation, triggering and additional cameras can be omitted.However, this relies on a periodically random behaviour of the water, which is normally the case.The methods can of course be applied on simultaneously acquired multi-camera datasets, too.It is not easy to favour one method over the other, because both have their pro's and con's.The object-space method requires a successful identification and measurement of underwater features in the single images.This would fail for images with heavily distorted image areas as shown in Figure 1.For such cases the image-space method seems to be without alternative.However, a tendency for loss of sharpness in the corrected image can be observed and therefore a loss of accuracy.In case of only moderate water surface movements, the object-space method should provide better results, as long as a feature can be identified and measured with a high level of accuracy in most images.Here a need for a modification of standard image matching methods can be identified, because the slight distortion of image features (due to refraction through wavy water) can't be described fully with a set of affine parameters.More investigations have to be carried on real world data at this point.

Figure 1 .
Figure 1.Disturbing imaging effects caused by glint and refraction.The grid pattern is displayed by a projector through the wavy water surface on to the water bottom.

Figure 2 .
Figure 2. Corrected image generated from an image sequence.The image content corresponds to an image of a projection through a calm planar water surface.

Figure 3 .
Figure 3. Experimental setup B (2 Nadir Cameras + 2 Oblique Cameras) with a tub filled with stones and water.The coded markers are for calibration purposes.

Figure 4 .
Figure 4.The measurement object imaged through calm water surface.Some of the river stones were marked with white crosses.Note the specular reflection from a ceiling lamp.

Figure 5 .
Figure 5. Oblique image of the moving water surface.Note the water pump (yellow circle) and the more crippled water surface nearby.

Figure 6 .
Figure 6.The range of movement of image features caused by a wavy water surface over a time period of 60sec.The length of the axes of the crosses marks the maximum deviation from the reference positions.The water pump is located in the middle of the lower edge of the scene.Therefore, the amplitude of feature movement is highest in this area.The underlaying image was taken through calm water.

Figure 7 .Figure 8 .
Figure 7.The Histogram of min/max movement of features over a time period of 60sec (60 images).

Figure 9 .
Figure 9.The Histogram of residuals [pixel] of positions of features in the corrected image vs. reference (image taken through calm water).

Figure 11 :
Figure 11: Simulated wave (Tessendorf model) and image data.The image coordinates of submerged points were calculated via ray tracing.The short yellow lines mark the surface normal in the point of intersection (yellow circle) of the image ray with the water surface.

Figure 12 .
Figure 12.Movement of an intersection point of an image ray with the water surface in Z over time (mean Z=0m, min/max = -0.4/0.5m).

Figure 13 .
Figure 13.Point Movement in image space (camera 1) over time due to wave effects.The red cross marks the point position correspondent to a calm planar surface.

Figure 14 -
Figure 14 -The ray tracing for one point and four cameras in two different water states.

Figure 15 -
Figure 15 -Bundles of image rays from 10 images of a sequence taken for evaluation of best fitting intersection point.The quadruple of rays, normally from different images, with the smallest distance from each other are marked in red.

Table 1
. RMS of calculated underwater point coordinates against the reference for different multi-ray setups