MITIGATING IMAGE RESIDUALS SYSTEMATIC PATTERNS IN UNDERWATER PHOTOGRAMMETRY

Systematic errors may result from the adoption of an incomplete functional model that is not able to properly incorporate all the effects involved in the image formation process. These errors very likely appear as systematic residual patterns in image observations and produce deformations of the photogrammetric model in object space. The Brown/Beyer model of self-calibration is often adopted in underwater photogrammetry, although it does not take into account the refraction introduced by the passage of the optical ray through different media, i.e. air and water. This reduces the potential accuracy of photogrammetry underwater. In this work, we investigate through simulations the depth-dependent systematic errors introduced by unmodelled refraction effects when both flat and dome ports are used. The importance of camera geometry to reduce the deformation in the object space is analyzed and mitigation measures to reduce the systematic patterns in image observations are investigated. It is shown how, for flat ports, the use of a stochastic approach, consisting in radial weighting of image observations, improves the accuracy in object space up to 50%. Iterative look-up table corrections are instead adopted to reduce the evident systematic residual patterns in the case of dome ports.


INTRODUCTION
Image residuals from bundle adjustment not only result from errors in the observations (random errors and gross errors), but also from errors in the functional or stochastical model (systematic errors) (Vlcek, 1969). When recognizable and repeatable residual patterns are identifiable in the image plane, this bias can be ascribed to residual systematic errors not properly modelled by the adopted mathematical model. The standard systematic errors that go beyond the collinearity conditions are well understood. The Brown/Beyer (Gruen & Beyer, 2001) model of self-calibration compensates those errors of amateur cameras and the Gruen model (Gruen, 1980) for professional photogrammetric cameras in in-air applications quite well. However, in underwater applications, because of the multi-media situation and additional optical system components, the systematic pattern has additional elements which are not covered by the conventional self-calibration parameters. The residual systematic errors are due to the incomplete functional model adopted, which alone is not sufficient to model the different refractive effects originating at the interface between water and the glass from a number of different sources such as, for example, a non-centered dome port, the use of a Ivanoff-Rebikoff corrector, a wet lens converter, or a flat port. This calls either for (a) additional self-calibration parameters or (b) a technique which was used in photogrammetry before self-calibration became a standard, called "Masson d'Autume method" (d'Autume, 1972) of post-adjustment systematic error treatment (see also Schilcher, 1980). Because of its simplicity in basic concept and ease of software implementation the method (b) was used here. Systematic image residuals are of concern since they may cause systematic errors in object space, often in the form of 3D polynomial deformations (Grün, 1978, James & Robson, 2014 * Corresponding author Nocerino et al., 2014). A compensation of these errors is only possible if they show up as residuals in either image space or at GCP residuals in object space. Otherwise they will deform the object space in a way that is non-detectable. For the mathematical proof see Grün (1986). When the standard self-calibration model is not sufficient, other approaches, such as collocation (Rampal, 1976) or finite elements (Lichti & Chapman, 1997;Tecklenburg et al., 2001) may be implemented. These approaches are computationally expensive and only little experience is available. On the other hand, an incomplete functional model can lead to systematic residuals in the image plane that are not only dependent on the position of the image observations within the sensor but are also a function of the 3D positions of the points in the object space. This is for example the case of image residuals for images acquired using a moving camera equipped with a rolling shutter sensor (Vautherin et al., 2016;Geyer et al., 2005) or submerged points in multimedia photogrammetry (Kotowski, 1988). In these cases, the residual patterns on the image sensor after bundle adjustment do not show clear systematic effects across all the images, and the residual function itself is dependent on the distance between the object point and the camera's perspective center (Maas, 1992;Li et al., 1997;Telem andFilin, 2010, Mulsow, 2010). In a previous study , the authors presented a stochastic approach that improves the accuracy of 3D reconstruction when using standard photogrammetric formulation underwater. A weighting function that penalizes the image observations more affected by optical aberrations and departure from the pinhole camera model has been used in a standard bundle adjustment computation, providing accuracy improvements up to 50% in real cases. In Nocerino et al. (2019) the authors showed that in the context of monitoring coral reefs through underwater photogrammetry performed by SCUBA divers, digital cameras mounted in dome ports display residual systematic patterns after bundle adjustment, which could not be compensated by standard selfcalibration functions. The residual patterns are ascribed to potential dome port misalignments. Through a redundant and geometrically strong camera network including nadir and oblique images, sub-centimetric accuracy results were verified against geodetic underwater measurements. On the other hand, in Piazza et al. (2018), in the context of photogrammetric monitoring of benthic species in Antarctica by SCUBA divers, the authors show how systematic errors remain uncompensated when using simplified protocols of image acquisition, by significantly bending the photogrammetric strip acquired along a simple transect. In this paper a novel analysis is presented aimed at quantifying and mitigating the effect of uncompensated systematic errors caused by flat ports and misaligned dome ports. In underwater photogrammetry, in particular by SCUBA divers, understanding, quantifying and controlling the errors within a limited budget, according to project tolerances, is fundamental, especially considering the complexity of the operations underwater where ground control is rarely available (limited time, high costs of operations, low accuracy). This can be accomplished only by an investigation where all the possible error sources are isolated and their contribution to the total budget assessed. With different camera housings, port types, size and adapters available on the market, it can be puzzling for nonexperts to understand whether a specific equipment combination is suited for photogrammetric analyses. Moreover, in long term ecological research projects, spanning over a period of several years, different photographic equipment can be utilized and interchanged (i.e. new cameras and housings or equipment failure); thus, knowing the potential systematic effects caused by possible misalignments of the dome port or use of a port with respect to another has very practical importance. With respect to the previous studies , where the combined effect of image quality degradation underwater and geometric errors caused by refraction were investigated in real cases, hereafter the two error sources are split and we focus only on the refractive effects through simulations. Simulations have the advantage that reference data is available for accuracy checking, which is otherwise very hard to generate underwater and that individual error components can be separated easily. The simulations are built upon a real underwater photogrammetry dataset collected within the Moorea Island Digital Ecosystem Avatar (IDEA) project (https://mooreaidea.ethz.ch/). A Nikon D750 full frame camera equipped with a 24 mm is used to run simulations with the same camera in a flat and misaligned dome ports where the entrance pupil of the lens does not coincide with the dome center. Errors in object space are then analyzed and mitigation strategies are proposed. These comprise different camera network image acquisitions and two photogrammetric processing strategies implemented by the authors in DBAT software (Börlin & Grussenmeyer, 2013) based respectively on radial weighting  and the use of an implicit method (d'Autume, 1972) integrated in an iterative selfcalibration approach for the correction of systematic image residuals. Experiments are reported in the case of self-calibration using standard underwater photogrammetry collinearity model with additional parameters.

MODELLING THE WATER EFFECT
Underwater, refractive effects can introduce systematic errors that are not considered by the standard pinhole model due to a number of factors that may include, to cite a few 1) the use of a flat port; 2) the non-perpendicularity of the flat port with the optical axis of the camera; 3) the non-concentricity of the entrance pupil of the lens with the center of the dome port; 4) the non-sphericity of the dome port. In this section the effects of systematic errors introduced using the simple Brown/Beyer functional model on accuracy potential underwater are investigated through simulations both in the case of flat and dome ports. A pinhole camera with principal distance and no distortions is placed in a pressure housing mounting, alternatively a flat and a dome port.

Modelling the refraction for a flat port
Referring to Figure 1-a, assuming that the flat port is at a distance ′ ̅̅̅̅̅ from the entrance pupil of the lens and that the optical axis is orthogonal to the flat port, the following equations hold: The International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences, Volume XLIII-B2-2020, 2020 XXIV ISPRS Congress (2020 edition) Due to the presence of water, the submerged point P is projected on the sensor at the distance ̅ from the principal point following the blue path in Figure 1-a according to Snell's law (equation 3).
If the camera were not immersed in water, the red collinearity line would instead directly link the object point with its image projection on the sensor, differing by the quantity ∆ ̅̅̅ with respect to the submersed case. The error ∆ ̅̅̅ changes depending on the position of point along the segment ̅̅̅̅ . The closer the point to the camera, the larger the ∆ ̅̅̅ variations. On the other hand, for small ′ ̅̅̅̅̅ compared to ̅ , the entrance pupil can be considered on the flat port, thus it may be assumed that the points and ′ coincide. In this case the ∆ ̅̅̅ is only function of the angle and it can be demonstrated that the error ∆ ̅̅̅ is very well modelled by a standard Brown/Beyer formulation .

Modelling the refraction for a dome port
As shown in previous works (She et al., 2019;Menna et al., 2016;Kotowsky, 1988), when the lens entrance pupil coincides with the center of the spherical dome, the refractive effect of water is negligible and collinearity is maintained. Contrarily, assuming that the center of the dome port with radius ̅̅̅̅ and the entrance pupil are misaligned by a distance ̅̅̅̅ along the optical axis, the following equations hold: From Figure 1-b it is clear that an axial misalignment with the center of the dome shifted towards the object space increases the actual field of view of the camera. On the contrary, the field of view is reduced when the dome center is behind the entrance pupil towards the camera sensor. Similar to 2.1, collinear points on the line ̅̅̅̅ produce different errors ∆ ̅̅̅ on the image sensor. The closer the point to the camera, the larger the ∆ ̅̅̅ variations.

Simulations
A real underwater photogrammetric image dataset was selected among those acquired in the context of coral reef monitoring long-term ecological research in Moorea, French Polynesia. This dataset is selected to build simulations with a camera network geometry that exemplifies real case acquisition, allowing to analyze the systematic errors introduced by refractive effects of both flat and dome ports. The images were acquired by a SCUBA diver at a depth of about 12 m with a Nikon D750 24Mpx full frame DSLR camera (pixel size ca 6 µm) mounting a 24 mm prime lens in a NiMAR pressure housing (Figure 2-a). The surveyed seabed area consists of three adjacent coral reef plots measuring 5x5m 2 each for a total area at the seabed of about 20x10m 2 , including a slightly larger buffer area around them. The surface characteristic is rough and hilly, with depth variation within the surveyed area of about 1.6 m (Figure 2-b). The relative camera to seabed distance was about 1.5 m with an average GSD of 0.4 mm. The dataset consists of approximately 600 images, of which 495 are nadir images, taken along 8 strips and one perimetral loop strip of oblique images pointing inside the plots (Figure 2-c). The dataset was processed using standard photogrammetric procedures, providing sub-centimetric accuracy, validated by using underwater geodetic surveying reference measurements (Nocerino et al., 2019). From the photogrammetric processed data, a subsampled 3D point cloud (1 point/20cm) was extracted and used as input to run the simulations. Three image networks were analyzed: 1. Complete network (nadir + oblique, 600 images). 2. Nadir network (500 images). 3. Single round-trip nadir strip (130 images).
The bundle adjustment with self-calibration on this dataset provides reference figures for the expected potential accuracy using the three different networks (complete, nadir, round-trip nadir). Using equations (1-13), image observations were generated to simulate the three datasets with misaligned dome and flat ports. A refraction index of 1.34 was chosen for water. The port thickness was neglected and the lens distortions were fixed to zero. In order to analyze the sole effect of the systematic errors induced by the misaligned ports, white noise was not added.
To summarize, four simulated datasets were created: A. Pinhole24, pinhole camera, no ports, no water. B. FP30, flat port at 30 mm from the entrance pupil. C. DP+30, misaligned dome port with dome center 30 mm ahead the entrance pupil . D. DP-30, misaligned dome port with dome center 30 mm behind the entrance pupil (towards the sensor).
The chosen offset values of 30 mm both for flat and dome ports are quite large. Although still possible, for example in the case of using wrong extension tubes, they are rarely seen in real underwater photographic equipment. In this simulation they were chosen as limit values to stress the related systematic errors. Starting from initial approximations and using the simulated image observations, minimal constraints bundle adjustments (BA) with self-calibration were performed in DBAT. As the systematic errors on the simulated dataset are expected to have only a radial component in image space, according to equations (1-13), bundle adjustment solutions were all computed only with radial distortion additional parameters. As accuracy assessment, errors were computed as the difference between the coordinates of triangulated image observations (sparse point cloud of tie points) and the input 3D point cloud through a similarity transformation. Table 1 reports the results of the self-calibration for the dataset A. From these results, it can be seen that, assuming a 0.25 pixel marking precision in image space, the focal length is estimated with a precision of respectively 0.25 µm (1/24 pixel) for the complete network, 1 µm for the nadir ones (1/6 pixel). Distortion coefficients are not statistically significant and the RMS of residuals on the image do not change across the three different networks. Figure 3 reports the result of precision in Z (the weakest measurement direction) for each tie point. The implemented image acquisition protocol with nadir and oblique images used for the coral reef monitoring in the Moorea IDEA project, is well suited for millimetric monitoring purposes with the entire plot, well below 0.5 mm. On the other hand, the single round-trip nadir strip results as the least accurate. Table 2 reports the results of the self-calibrations for the datasets B, C, D. In these cases, the RMS of image residuals are only due to the incomplete functional model, unable to properly take into account the refraction effects of water. The RMS of image residuals for the flat port are significantly larger (almost ten times). The nadir camera networks provide three times smaller RMS residuals in image space. Although the RMS of image residuals is often considered a parameter of internal precision, low values do not correspond to higher accuracy per se. Indeed, as it can be seen from Figure 4, weak camera network geometry tends to absorb systematic errors within the exterior orientation parameters due to projective couplings, with camera positions that move and rotate to reduce the image residuals (see also Grün, 1986). Indeed, looking at Figure 4, although the single round trip strip records the best RMS for image residuals, it is also the one with worst accuracy in object space displaying the largest RMSEZ among the same dataset. Self-calibrations for the flat port are always less accurate of almost an order of magnitude with respect to dome ports. For the datasets flat port FP30 and dome port DP-30 (dome center between the entrance pupil and the sensor) the refraction introduces pincushion distortions with the maximum values for the flat port almost five times larger than for the dome port. Also, an increase of principal distance is recorded (field of view reduction); on the other hand, the DP+30 dataset shows a reduction of principal distance (larger field of view) and barrel distortion. A principal distance variation between calibrations inside and outside the water is thus a method that can be used to fine tune the dome port centering. Figure 4 shows the color-coded error maps reporting, for each tie point of datasets B, C, D, the error in Z, the largest among the three components. Dome shaped deformation is visible for all datasets, a clear sign of uncompensated systematic errors. Complete camera networks confirm the results from previous studies (James & Robson, 2014;Nocerino et al., 2014) with significant mitigation measures and improvements of up to three times. 12.0002 ± 5.6e-06 12.0001 ± 3.08e-06 12 ± 6.88e-06 k1 [mm -2 ] 9.89067e-05 ± 2.6e-09 9.89603e-05 ± 1.52e-09 9.9078e-05 ± 2.99e-09 k2 [mm -4 ] -3.15199e-08 ± 1.25e-11 -3.04765e-08 ± 7.01e-12 -3.05642e-08 ± 1.36e-11 k3 [mm -6 ] 1.56546e-11 ± 1.83e-14 1.39155e-11 ± 1.  Table 2. Results of the bundle adjustment with self-calibration for the simulated dataset B,C,D.
As already highlighted in (Menna et al, 2017), the flat port shows the worst accuracy with a slight improvement when more redundant and complete camera networks are used. Going from the single round-trip strip to complete network, an improvement of 25% is obtained for the RMSEZ (from 8 mm to 6 mm) and of 36% for the maximum error span that reduces from 46 mm ([-14.6 mm, +31.7 mm]) to 30 mm ([-10.3 mm, +19.6 mm]). For the datasets C, D, despite a large dome centerentrance pupil axial misalignment of 30 mm, the refraction effect is very well mitigated by a camera network with both nadir and oblique images, providing a RMSEZ well below the millimeter.

MITIGATION OF RESIDUAL SYSTEMATIC PATTERNS
This section describes the implementation of computational methods for mitigating unmodeled systematic errors in a selfcalibrating bundle adjustment approach. Both methods have been implemented by the authors in MATLAB, currently in the form of a plugin for DBAT. The first method  uses radial weighting (RW) of image observations which proportionally penalizes those observations with larger radial distances. The second method estimates look up table corrections from systematic residual patterns iteratively. After the first selfcalibrating bundle adjustment solution, the image plane is subdivided into a grid, where for each cell the median image residual is computed for x and y image coordinates. The correction, equal to the computed median residual error, is applied to the x and y image observations and a new selfcalibrating bundle adjustment step is run. The procedure is iteratively repeated until the convergence criterium is reached (difference in the solution vector).

Radial weighting
The dataset B FP30 was processed using a radial penalty implemented through a weighting function in DBAT . A standard deviation of 0.01 pixel was given to image observations at the center of the image format growing up to 30 pixels for the image corners following a third-degree polynomial. Table 3 shows significantly different calibration parameters and a slight improvement of their precisions. Figure 5 shows that the accuracy in the object space improves up to 50% despite an increase of the RMS of image residuals; this behavior confirms previous findings obtained in real case studies underwater. For the dome port datasets, improvements using radial weights were not significant as for the flat port and results are not reported in this contribution.

Residual systematic pattern correction in simulations
The technique we implement here uses two distinct look up table corrections using the median values computed within square cells whose size is typically chosen between 64x64 and 128x128 pixels 2 depending on the pattern behavior (spatial frequency). In the current implementation the cell size is chosen manually.
To demonstrate the working principle of the implemented technique, we first show how it works for a case where the residual systematic pattern is very evident. The single strip dome port dataset is processed with only k1 as unique additional parameter to absorb the refractive effects of water. Figure 6-a shows the classical representation of residual vectors, here computed as average within a grid with cells of 64x64 pixels 2 ; on the right column the image residuals are plotted against their radial distance; positive residual magnitudes mean that the residual vectors point towards the corner and negative ones towards the center of the sensor. Although residuals of opposite directions exist at each radial distance, a systematic harmonic trend can be seen in both representations. Figure 6-b shows the results after applying the implemented procedure. The magnitude of residuals is reduced (RMS lowered from 0.372 to 0.029 pixels) and a systematic component is removed with an improvement of RMSEZ from 0.025 m to 0.0011 m, obtaining values similar to when the complete set of additional parameters k1, k2, k3 are used.   Table 3. Results of the bundle adjustment with self-calibration for the simulated dataset B, using radial weighting (RW).
After applying the proposed method, the systematic harmonic pattern (Figure 6-a) is flattened and reduced in magnitude ( Figure  6-b). Nevertheless, it is worth noticing that Figure 6-b left reports the average residuals computed within the grid cells: residual values equal to zero at specific points on the sensor do not necessarily indicate that the residuals have been completely corrected, but rather that their average has been zeroed. Indeed, from Figure 6-b right, for a specific radial distance, both positive and negative residuals are visible, although their average has become zero after the correction, which has effectively removed the systematic component. The method is also applied to the dome ports with the complete set of k1, k2, k3 radial distortion parameters (Figure 6-c). In this case, the systematic effects are removed as seen in Figure 6-d left after only one iteration, with a reduction of the RMS of image residuals from 0.009 to 0.008. Nevertheless, no significant accuracy improvements are found in object space. This behavior can be ascribed to the fact that the magnitude of systematic pattern residuals is very small (hundredths of a pixel) and has no practical effects on improving 3D point accuracy in object space for the analyzed dataset.

DISCUSSION AND CONCLUSIONS
The presented simulations showed the effects of systematic errors caused by unmodelled refraction and use of a standard photogrammetric functional model. With respect to theoretical expectations for the potential accuracy without modeling the refractive effects, a worsening of up to a factor 40 is observed for the flat port and a factor 2 for the dome port. The use of a redundant camera network with both nadir and oblique images significantly mitigates the systematic errors, with complete networks being always better than just nadir ones by a factor of at least two. For the flat port, the use of a stochastic approach, radial reweighting of image observations, confirms an improvement of accuracy up to 50% in object space. The refractive effects cause clear systematic residual patterns in the case of dome ports. For the dome port, the use of mitigation measures based on lookup table corrections, iteratively computed within a bundle adjustment approach, did not significantly improve the accuracy in object space due to the low absolute values of the corrections.
The International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences, Volume XLIII-B2-2020, 2020 XXIV ISPRS Congress (2020 edition) a) b) c) d) Figure 6. Image residuals for the dataset DP+30 SINGLE ROUND TRIP STRIP before (a,c) and after (b,d) applying the developed method. In (a) only k1 is computed in the bundle adjustment; (b) shows the results after further applying 10 iterations to (a) with the proposed method. Full radial distortion parameters k1, k2, k3 are computed in the bundle adjustment (c); (d) shows the results after further applying 1 iteration to (c) with the proposed method. Further developments of the implemented approaches, currently under investigation, consist of the combination of the two methods and the use of a four dimensional look up table where the third input parameter is the relative depth of object points. Figures 7a and 7b show image residuals versus the relative depth for a grid cell of 64 pixels at two specified radial distances of 1184 and 2784 pixels for the single round trip strip dataset DP+30. A clear depth dependent pattern is observed.