Multimedia Photogrammetry with non-planar Water Surfaces – Accuracy Analysis on Simulation Basis

If multimedia-photogrammetry is used for the generation of point clouds of submerged objects or of the water bottom, Snell’s law has to be considered. When the images are taken from air, image rays are refracted at the air-water interface. This results in the collinearity equations being no longer valid. Bundle block adjustment can still be solved by adding additional terms considering Snell’s law. Existing approaches usually assume that the water surface is flat. Refractive indices and water height can either be measured separately or included as unknowns in the adjustment. However, when the water surface is not flat due to the presence of waves, assuming a planar water surface leads to large geometric errors. This work will analyze the significance of those errors and propose a way of including water surface parameters as unknowns into the bundle block adjustment, both based on simulated data. The simulation reproduces multiple images taken simultaneously, e.g. from synchronized UAV cameras or from cameras on tripods.


Introduction
Multimedia photogrammetry has been a topic of research for the last years and decades (Maas, 1995;Mandlburger, 2022).It can be divided in two groups.In the first, the camera is carried underwater, i.e. by a scuba diver or an underwater vehicle.In this case, the camera is placed in a waterproof housing with either a flat or a dome interface, compensating the influence of refraction on the air-glass-water interface.Using a (hemi-) spherical interface enables a geometry, where all image rays pass through the interface perpendicularly.In this case, the image rays are not refracted (Menna et al., 2016).Using a coplanar flat interface reduces the refraction influence to an effect that can largely be absorbed using conventional sensor modelling and camera calibration methods (Shortis, 2015).Placing the camera oblique to a (planar) interface requires strict geometric considerations of Snell's law (Mulsow, 2010).Underwater image data acquisition is for instance applied in archaeology (Menna et al., 2018) or 3D measurement of aquatic ecosystems (Remmers et al., 2024).The second group of multimedia photogrammetry applications uses cameras that are placed outside (above) the water body.Thus, it is also denoted as through-water photogrammetry.Herein, cameras are for example placed on a tripod or on an uncrewed aerial vehicle (UAV).Refraction occurs on the water surface directly from air to water.Due to the difficulty of measuring wave water surfaces, in many applications, the water surface is assumed flat and horizontal.For the determination of the geometry of submerged objects or the ground, Snell's law has to be considered separately for every image ray that travels from the image sensor through the camera's projection centre and the air-water interface to the object point.For the consideration of a flat, horizontal water surface, this is reduced to two parameters: water height and refractive index.Potential applications include laboratory experiments or bathymetric measurements from UAVs (Mulsow et al., 2020).The approach of this study addresses a problem that occurs in through-water photogrammetry, where the interface between air and water is not a stable glass housing but the open water surface, which will be non-planar in most cases.In case of wind, currents, or other (possibly also artificial) influences, the water surface becomes wavy or in other ways non-planar.This leads to refraction influences being more complex.Image rays are refracted individually, depending on where they intersect the water surface.The intersection points vary in surface height and surface inclination.Every image observation provides two measurements (image coordinates x' and y').An object point has three parameters (X, Y and Z) that need to be determined.If camera orientations (interior and exterior) are assumed to be known and an object point is seen by n cameras, there are 2n observations and 3 parameters.This leaves 2n-3 degrees of freedom for the water surface.If we try to determine the intersection points for every image ray, in order to calculate the refraction using Snell's law, we need three position and two orientation parameters for every intersection point.This results in a negative degree of freedom of -3n-3 and the problem is therefore underdetermined.Thus, it is necessary to model the water surface with an appropriate number of parameters that ensures determinability on the one hand and describes the water surface in sufficient detail that the refracted underwater image rays intersect in the submerged object point on the other hand.

Simulation Process
To evaluate the influence of non-planar water surfaces in through-water photogrammetry, a simulation pipeline was established.All steps are performed in MatLab using a combination of existing toolboxes plus our own functions.The simulation includes three parts: The first part contains the generation of synthetic images using pixel-wise forward raytracing.The second part contains the projection of specific object points into image space.The third part contains the bundle block adjustment.The three parts are described in more detail in the following paragraphs.

Rendering of Synthetic Images
Rendering a synthetic image needs three main input pieces of information: (exterior and interior) orientation of the image, a representation of the water surface, and a representation of the underwater object space.The water surface is represented either globally as a functional equation or piecewise locally using splines, both being described in Section 3. The object space may be represented as a point cloud with intensity information (i.e.RGB) for each point.To render a synthetic image, the 3D image ray in air  ⃗  is calculated for every pixel of the image using the interior and exterior orientation parameters of the image.Subsequently, the intersection point with the water surface has to be found.In most cases, this is an iterative process for both water surface representations (see Section 3).Together with the surface normal N at that intersection, a 3D vector form of Snell's law is applied to receive the underwater image ray vector  ⃗  : Applying this method to all pixels of an image, enables the rendering of a synthetic image.The image quality and the computing time both depend on the number of pixels (resolution) and the density of the object point cloud.
Figure 1.Refraction of a 3D image ray  ⃗ ⃗⃗  (black) that intersects the water surface with the normal  ⃗⃗⃗ (red).The vector is refracted according to Snell's law and continues under water as  ⃗ ⃗⃗  (green)

Projection of Object Points into Image Space
Projecting object points into image space with known exterior and interior orientation of the camera is a core task of photogrammetry, usually directly solved with the collinearity equation.In the multimedia case, the image rays have to be tracked from the object point through the interface to the image space including Snell's law.Mulsow et al. (2010) introduced the term of alternating forward ray tracing (AFRT) as an iterative realisation of backward ray tracing.In this simulation, we use AFRT.For a first approximation, the closest integer pixel position from the synthetic image is used.Alternatively, the approximate image point could be calculated with the collinearity equation, assuming no water being present.Starting with that approximation, forward ray tracing through the water surface is applied and the 3D distance from the underwater image vector to the object point is calculated.The x' and y' coordinates are then iteratively optimized in sub-pixel-space until that distance reaches a minimum.

Bundle Block Adjustment
For bundle block adjustment, we follow the approach of Rofallski and Luhmann (2022), where minimization is conducted in object space rather than in image space.Finding homologous points in multiple images could be done with feature point detectors and descriptors on the synthetic images generated in 2.1.However, in order to avoid inaccuracies in image space and the handling of mismatches, we use error-free projections of given object points (2.2) in multiple images.There are three groups of parameters that can be determined in that adjustment: (interior and exterior) camera orientation parameters, object point coordinates and water surface parameters.For this study, we assume camera orientations to be known.Interior orientations can be calibrated in advance.Exterior orientations are either fix (in a laboratory environment) and can also be determined separately, or they may be determined from other sensors (GNSS and IMU) or from parts of the image that show above-water objects.This leaves three combinations of parameters to determine in the adjustment: object points, water surface, or both.The specific challenges of those combinations are described in the following paragraphs.

Known water surface and unknown object points
In the first version, we assume a known water surface, with one of the representations from Section 2. Using error-free image points and correspondences, object point coordinates can be determined in a spatial intersection including refractions at the air-water interface.This case does not reflect many real world application, but is still interesting to simulate to get a better understanding of the model.The functional model of the bundle block adjustment has to be extended with terms for refraction.We furthermore design the bundle block adjustment to minimize the distances between the underwater image rays and the object point (as in Rofallski and Luhmann, 2022).The underwater image ray has to be determined from the intersection point with the water surface and the surface normal at that point.These are calculated depending on the water surface representation (cf. Section 3).Since the water surface is known in this scenario, it does not add any unknowns to the adjustment.However, it still does add computation time, since the determination of the intersection with the water has to be solved iteratively.

Known object points and unknown water surface
If the underwater object point coordinates are known, this version can be used to determine the water surface geometry.This is for example relevant in laboratory experiments as in Rupnik et al. (2015).The coordinates of underwater markers has to be known, for instance as a checkerboard pattern.They can for example be measured before water is added.Determining the water surface in the adjustment adds parameters to the adjustment.The number of parameters depends on the utilized water surface model.It has to be considered that the system of equations remains overdetermined.The available degree of freedom depends on the number of images and homologous points.The number of parameters in functional equations or knots in splines therefore needs to be adapted accordingly.Furthermore, approximate values for the parameters are required.The model is especially sensitive to the choice of approximate values if the water surface is defined with a global function.

Unknown water surface and object points coordinates
Combining the determination of water surface and object point coordinates in a common adjustment is the most sophisticated of the presented methods.The redundancy is further reduced, necessitating either a small number of parameters to describe the water surface or a large number of well-distributed image observations.This study only covers the determination of either water surface or water bottom.This forms the foundation for the simultaneous adjustment of both, which will be object of future research.

Water Surface Representations
In this study, we use two methods of representing the water surface: with a global functional equation and with cubic splines.Both methods as well as the determination of the intersection point with an image ray are described in the following paragraphs.

Functional Equations
The first water surface representation uses a functional equation: where: Xw, Yw, Zw = water surface coordinates P = water surface parameters A simple example for this is a sine wave along the X axis, which is defined as: This can be solved using the Newton approach.The normal at that position can be determined from the derivatives: where:  = derivative in x  = derivative in y

Cubic Splines
The second water surface representation can be used for a wider range of water surface geometries.A freeform geometry is hereby interpolated by cubic splines.The splines interpolate a given cloud of knot points.For the determination of the intersection between image ray and water surface, it is important that the interpolation function is universally differentiable, which is the case for cubic splines that are C2 continuous.The intersection of a vector with the splines is solved in an iterative adjustment, minimizing the distance of the intersection point candidate to the splines and to the image vector.The normal is calculated with Equation (5).For the determination of an unknown water surface, a regular raster of knots can be used that unambiguously defines the cubic spline interpolation.Their Xand Y-coordinates are kept fix, while the Z-coordinates are determined in the adjustment.

Simulated Environment
The simulation for this study consists of an underwater checkerboard and four nadir looking cameras with approx.60% overlap in X and Y.The size of the checkerboard and the position of the cameras are given in Table 1.The simulation is designed in a way that the median water depth is 1 (unit-free) and can be scaled arbitrarily, as long as the ratios stay the same.1. Simulation parameters For this study, two water surfaces were generated and imaged with the same camera orientations and checkerboard object.The first is a sinusoidal wave with a wavelength of 1.5 and an amplitude of 0.25.The second is a freeform water surface that was generated in the software Blender.A simulation of water flowing into a water basin was conducted in Blender and one frame of that procedure was used for the simulation.The freeform water surface was sampled with 7156 points (point spacing 0.1 in X and Y) and a cubic spline interpolation was fitted to those points.

Results
First, synthetic images for all four camera positions were generated using pixel-wise forward ray tracing as described in Section 2.1.The images are used for a better understanding of the influence of the non-planar water surface and to produce approximations for the sub-pixel image observations.Figure 3 shows the images generated from the blue camera position in Figure 2. In the second step, all crossing points of the checkerboard pattern were projected into image space and used as control points for further analysis.For simplification, the bundle block adjustment uses only those corner points that are visible in all four images.

Influence of Water Surface on Object Coordinates
When a commercial photogrammetric or structure-from-motion software is applied in environmental mapping using UAV images, underwater areas are treated just as all other areas, applying the collinearity equations.Mulsow et al. (2020) showed a workflow that adds the refraction influences that occur on a planar water surface.This is still insufficient, when the water surface is non-planar.Figure 4 shows the effect of assuming either no water or a planar water surface on the determination of an object point from two images.The sinusoidal and freeform water surfaces introduced in Section 4 were used to simulate the error that occurs under the false assumptions of no water or a planar water surface.Therefore, the checkerboard corners where projected in the four images using the approach described in Section 2.2 and the according water surfaces (sine wave and spline).This simulates images taken through a non-planar water surface.The image observations were then used for object point determination assuming no water, a planar water surface (at mean water height) and the correct water surface.The calculated object points were compared to the ground truth object point coordinates and a root mean squared (RMS) 3D-distance was calculated: where:  0 = ground truth object point coordinates  = calculated object point coordinates ‖ 0 − ‖ = 3D-distance  = number of object points The results are shown in Table 2.The zero RMS for the correct water surface proofs that the intersection with a non-planar water surface (section 2.3) works correctly.The other values underline that assuming a planar or no water surface leads to considerable errors.
The effect of falsely assuming a planar water surface depends on the shape of the actual water surface and on the depth of the object points.We therefore also tested sinusoidal water surfaces with different wavelengths and amplitudes in the simulation environment and calculated the RMS that occurs, when a planar water surface is assumed (Figure 5).Shorter wavelengths result in steeper waves.The surface normals therefore differ decisively from a planar water surface, producing larger errors, when a planar water surface is assumed instead of the correct wave.
Greater amplitudes also result in steeper waves.The intersections of the image rays with the water surface vary more in height, also resulting in an increasing RMS when a planar water surface on the mean water height is assumed.Additionally, we varied the water depth of the checkerboard pattern in the simulation for the initial sinusoidal water surface (Figure 6).The RMS increases with depth.
Figure 5. Influence of a sinusoidal water surface on object point accuracy depending on wave parameters wavelength (left) and amplitude (right), when a planar water surface is falsely assumed.The parameters that are not varied can be found in Table 1.
Figure 6.Influence of a sinusoidal water surface on object point accuracy depending on water depth of checkerboard points, when a planar water surface is falsely assumed.The parameters of the water surface can be found in Table1.
For freeform water surfaces, the accuracy depends on the resolution of the spline interpolation.The water surface was initially interpolated with cubic splines based on approx.7200 points.In order to analyse the influence of the interpolation resolution, we sampled a varying number of points (regularly distributed) from the originally modelled freeform water surface.By interpolating the sample points with cubic splines and integrating them into the object point determination, the accuracy in proportion to the spline resolution can be determined (Figure 7).It can be observed that the accuracy depends on the density of the sample points.Using 4 knots reproduces a result that is only slightly better than the planar surface (Table 2).With an increasing number of knots, the accuracy first increases almost linearly.It can still be observed, that even a large number of regularly rasterized knots is not able to reproduce the exact result.
Figure 7. Influence of spline resolution on object point accuracy.The freeform water surface from Table 1 was sampled with regular distributed knots, interpolated by cubic splines and used for object point determination.

Determination of Water Surface from Control Points
When object point coordinates and image orientations are known, the water surface parameters can be determined in an adjustment on the basis of water bottom point image coordinate measurements.This was tested with the sinusoidal and the freeform water surface.Both water surface parameterizations have their own challenges.For the sinusoidal water surface, there are two main challenges: The first one is caused by the global definition of the water surface.A small change in wavelength for example adds up to very large differences in the appearance of the wave further away from the defined origin of the wave.This demands the adjustment to start with rather good approximate values.The second challenge lies in the repetitive behaviour of the wave.The search for the correct intersection of an image ray with AFRT or backwards ray tracing tends to alternate or diverge, as already described in Mulsow (2010).The challenge of the freeform water surface lies in the resolution of the interpolation.As shown in Figure 7, the accuracy depends on the number of points that were used to sample the water surface.

Summary and Outlook
This study presents a simulation environment that addresses the topic of through-water-photogrammetry with non-planar water surfaces.It was used for the generation of synthetic overlapping images of a checkerboard pattern recorded through a sinusoidal and a freeform water surface.Projections of the checkerboard corners in four images were calculated using backwards ray tracing (realised by alternating forward ray tracing).The synthetic image observations were used to simulate the influence of the water surface geometry on the object point determination.
We showed that assuming a planar water surface leads to considerable errors depending on the shape of the water surface and the imaging geometry.When object point coordinates are known, it is possible to determine the water surface geometry in a bundle block adjustment, with requirements for the approximate values (functional) and limitations in resolution (freeform).Future research will address the simultaneous determination of water surface and underwater object coordinates.This will be done both in the simulation environment and on real data.
Figure 9 shows a laboratory experiment with four cameras recording a water tank where the water is brought to movement by a small water pump, generating sine-like waves.
Figure 9. Lab experiment with four synchronized nadir cameras that observe a tank with water that is put in motion by a water pump.
⃗ ⃗⃗ ⋅  ⃗  ) − 1))  ⃗ ⃗⃗ (1) where  ⃗ {,} = ray vectors in air and water  {,} = refractive indices of air and water  ⃗ ⃗⃗ = normal of interface at intersection This vector  ⃗  is then intersected with the point cloud.The RGB values of the three closest points to  ⃗  are interpolated according to their distance to the vector and saved for the according pixel.
) where:  ̅ = mean water height  ̂ = amplitude  = wavelength The water surface can therefore be described with three parameters.The intersection point of an image ray  ⃗  = [   ,    ,    ] can be determined by inserting the image ray into the equation of the water surface and determining the zero point:

Figure 2 .
Figure 2. Simulation environment with freeform (top) and sinusoidal (bottom) water surface.The purple, yellow, red and blue lines show the intersection of the image borders with the water surface and the underwater checkerboard.

Figure 3 .
Figure 3. Synthetic images with sinusoidal (top) and freeform (bottom) water surface.The images were generated from the blue camera positions in Figure 2 (top and bottom).

Figure 4 .
Figure 4. Influence of water surface on object point position in 2D.The image rays in air (solid black) from two cameras (black circles) are refracted assuming a sinusoidal (blue) or a planar water surface (orange) resulting in different intersections under water (blue and orange x).The dashed black lines show the straight continuations of the image rays with no refraction.
For a good accuracy, a large number of sample points should be determined.The total number is on the other hand limited by the tolerable degree of freedom that depends on the number of image measurements.The splines interpolate the water surface locally.Each knot only influences a certain area of the water surface.The positions of sample points in areas where only a few or no image rays penetrate the water are therefore difficult or impossible to determine.We tested the determination of the freeform water surface from Figure2(top) with 256 knots, even when the number of image observations would have allowed for a denser interpolation.The Xand Y-coordinates of the knots are distributed regularly and kept fix.The Z-coordinates are approximated with the mean water height and determined using the Gauss-Helmert-model.The adjustment starts to alternate after 12 iterations.The distances from the resulting underwater image rays to the object points have an RMS of 0.025 at that stage of the adjustment.This is even superior to the simulation in the last paragraph of Section 5.1.Observing the resulting water surface reveals that the resulting water surface is not consistent with the ground truth water surface (Figure8top).The RMS of the adjusted Z-values on the 256 knots compared to the ground truth is 0.5.Large deviations mainly occur at the border of the interest area, which is covered only by a few image rays (Figure8bottom).The small RMS in object space suggests that the resulting surface might be the best mathematical solution for that input data, but does not reproduce the correct water surface in the edge areas.However, the centre of the water surface on the other hand has been reproduced much better.The central area (−1.5 >  > 1.5;−2 >  > 2) contains 36 of the 256 knots.The RMS of the adjusted Z-coordinates of the water surface knots in that area is 0.03.

Figure 8 .
Figure 8. Water surface (top) resulting from an adjustment of a cubic spline interpolation from known object points.The black points represent the crossing points of the checkerboard that are visible in all four images.The ground truth water surface can be found in Figure 2. The differences in Z on the 256 knots are shown in the bottom image.The central region is marked with a red rectangle.