AN EFFICIENT SOLUTION TO RAY TRACING PROBLEMS FOR HEMISPHERICAL REFRACTIVE INTERFACES

: Refraction effects, their description and modeling are important aspects of underwater and multimedia photogrammetry. For hemispherical interfaces, the usual approach to refraction is to rely on standard pinhole representations, e.g. by employing the Brown model. This is strictly only possible if entrance pupil of the lens and dome center coincide which is not trivial to achieve. However, simulations and other authors show that systematic residual errors occur with these approaches up to considerable margins if offsets of some millimeters are present. Hence, we propose a novel efficient, yet strict optimization algorithm to account for offsets between dome port centers and entrance pupil. It is about two orders of magnitude faster than standard ray tracing implementations that account for refraction while providing similar or equal results. The algorithm is employed for analysis on a simulation and two real data sets and performance of additionally estimating the dome center is investigated. Our method is capable of improving accuracy in one data set at a maximum of 30 % but even so cannot provide improvements for the second data sets. An explicit calibration model is hence to be chosen carefully and most likely relies on the offset’s margins and each individual application.


INTRODUCTION
Multimedia imagery, especially captured under water, suffers from many detrimental factors affecting both image quality and geometry, compared to single-media image acquisition. Firstly, the light from the object travels through multiple media (air, glass, water) and is refracted at the separating interfaces. This, among other influences, introduces astigmatism, chromatic aberration and absorbs colors from long wavelengths in shallow waters to medium and short wavelengths in higher depths. This creates color shift towards green and blue in medium depth waters to approx. 30 m and decreases image quality and sharpness (Shortis, 2015). Moreover, the ray path is altered at refracting interfaces according to Snell's law, thus rendering the standard pinhole model with additional distortion parameters invalid (Treibitz et al., 2012). For flat interfaces, this is often compensated by either calibrating under water with standard parameters of photogrammetric models (e.g. Brown 1971), designed for air-applications to implicitly calibrate the imaging geometry, or by explicit calibration of the refractive imaging properties.
Refraction can also be encountered by employing a hemispherical interface (dome port). If the entrance pupil of the lens is aligned with the hemisphere's center, all rays pass the interface orthogonally and no refraction occurs. However, it is a rather complicated task to mechanically align a camera inside a dome housing to the point that no residual errors are introduced (e.g. by visual servoing). Very fine positioning with special methods are required for that (She et al., 2019). Even when aligned, the position of the entrance pupil may change with varying focal distances which enforces realignment for different settings . This also applies to * Corresponding author added water pressure which might force dome ports out of their aligned position. As investigated by Kunz and Singh (2008), unresolved misalignments result in nonlinear distortions that cannot be fully compensated by standard distortion models. Another disadvantage is that the spherical shape of a dome port creates an inverse lens in front of the camera and hence virtually images objects at about three times the dome radius. This can produce very close projections with small dome radii of a few inches, as is the case with many standard underwater housings . If the camera lens is not capable of focusing at such short distance, a dome port cannot be employed to its full capabilities and blurry images are formed. This restriction could be omitted by adding a diopter for macro focus to the lens which would in turn include another element into the optical ray path Menna et al. 2017a).
One major advantage of an explicit dome calibration is, apart from a possible accuracy improvement, the possibility to align the dome port more accurately by applying the calibrated values to the dome position until no significant distortions occur. This could be performed in a closed control loop that could even be automated (e.g. visual servoing). When using non-customized dome ports, as e.g. from BlueRobotics 1 , it is rather complicated finding the dome center (Kahmen and Luhmann, 2022). With certain lenses, it can even become impossible to center the dome because the camera has to be moved very close to the dome which might be blocked by a sealing flange or the lens itself. Hence, it is desirable to be able to account for such cases by introducing an explicit calibration model. Dome ports are usually better suited for any kind of optical application when well-centered. This holds especially true in high pressure environments, e.g. deep sea, where the favorable spherical geometry is better suited, as forces spread wider on the entire surface (She et al., 2022). Additionally, the field of view is limited and refraction is significantly higher resulting in decreased image quality with flat ports (Menna et al., 2017a).
The most relevant contributions of this paper are an efficient optimization procedure for ray tracing with hemispherical dome interfaces on the one hand and the transfer of the procedure to analytical derivatives for a significant computation time improvement towards real-time applications on the other hand.

RELATED WORK
The alignment of a dome center with the camera's entrance pupil is not trivial because both the entrance pupil and the dome center are not physically tangible. For one, She et al. (2019) proposed exploiting the geometric image constraint that straight lines over and under water have to remain straight lines at the media transition if the alignment is correct. Hence, the authors set up a camera such that the dome is half-way under water and faces parallel to the water surface. The camera can be moved accurately in axial direction and hence iteratively positioned to collimate lines under and over water. This approach requires a specially constructed water container with an integrated dome and linear unit. This might be impractical for many operators to obtain. Kahmen and Luhmann (2022) proposed an individually designed additive shell construction that is mounted onto the front of a camera lens. The shell is an inverted shape of the dome and constrains the camera center to coincide with the dome center, given an accurate axial centering of the camera. The distance for the shell construction can be obtained from data sheets of the dome and the lens and have to be adapted for each individual lens-dome setup. Minor adjustments regarding focus settings may also be necessary since the entrance pupil shifts with focus, as stated by Menna et al. (2016). Results showed very high imaging quality but the authors have not quantified the accuracy of the alignment, yet. Kunz and Singh (2008) simulated the effect of both, a decentered dome port (on the order of one millimeter per direction) and a flat port with a simplified stereo baseline of 75 cm. The basic physics of Snell's law are applied iteratively to obtain image coordinates that fulfill the refractive geometry under water and a standard pinhole model is used to calculate 3D coordinates from the resulting image coordinates. The results show larger deviations on the order of centimeters for the flat port and smaller ones for a dome port of a few millimeters. However, both procedures produce systematic errors in object space with highest magnitudes in the image corners, if refraction is neglected. Contrary, Nocerino et al. (2016) provide an order of few millimeters where no significant distortions are said to occur. The same is stated for offsets from concentricity of the two dome radii.
References, such as Kunz and Singh (2008) or Nocerino et al. (2016) state that dome systems underwater can well be modeled by the standard pinhole model with additional distortion parameters, e.g. by Brown (1971). However, especially in highly accurate applications, such as underwater inspection (e.g. Kahmen and Luhmann 2022), calibrating dome offsets remains of interest to exploit the full potential of underwater photogrammetry. She et al. (2019) tackled the issue by taking an image in both air and water at the same position and optimizing for the residual errors with least squares. Later, in She et al. (2022), the authors extended the model of an axial camera, formulated in Agrawal et al. (2012) and Treibitz et al. (2012) for flat ports and applied it to the hemispherical dome case. With additional optimization of projected errors on a checkerboard, the need for air/water images of an object at the same position are not required anymore, creating a leaner calibration process without additional iterations. Iscar and Johnson-Roberson (2020) calibrated a stereo camera system inside a single dome port housing which creates major refractive effects in the imagery. The authors proposed, using the point spread function to describe the refractive effects. From there, correction on both geometry and the image itself can be performed. The procedure itself however, requires an extensive experiment of about 200 camera positions and constant environmental conditions which limits the approaches practical feasibility. Yang et al. (2021) developed a methodology, based on the coplanarity constraint, to account for spherical refraction. However, their application is based on mining environments with refraction only taking place at a rather small glass interface which is of considerable less effect than the underwater application with air-glass-water ray paths.
The ray tracing approach, introduced by Kotowski (1987Kotowski ( , 1988, and subsequent extensions based on the same approach, e.g. by Mulsow (2010) or Mulsow and Maas (2014), are generically formulated to cover all kinds of surface shapes of first and second order.
However, the works mostly focus on object-invariant cases with different shapes, such as cylinders and cuboids. Special considerations on dome calibration, such as the determinability of offsets from the dome center and refractive indices, are not covered. Rofallski and Luhmann (2022) improved the optimization strategy for the aforementioned ray tracing algorithms by shifting the error function from image space to object space. Hereby, computational speed was significantly increased while maintaining a strict physical representation of the imaging geometry for flat ports. Bosch et al. (2015) developed an omnidirectional camera and calibrated it, employing ray tracing through a sphere-cylinder composite. The procedure is iterative and comparable to Kotowski (1987) and Mulsow and Maas (2014) to reproject the object points to image space. However, explicit calibration can introduce high mathematical correlations between the interior orientation (i.e. camera constant c and principal point, x ′ p , y ′ p ) and the dome center coordinates. It may thus be advisable to calibrate in a two-step procedure, as recommended by Bosch et al. (2015).
This contribution bases on the findings by Rofallski and Luhmann (2022) for flat ports and extends the basic idea of optimization in object space for dome ports. Hence, we employ ray tracing fundamentals and include an improved optimization strategy to perform bundle adjustment for underwater photogrammetry with highly improved computation speed. Furthermore, we show that the model is capable of handling analytical derivatives which is a major shortcoming of other ray tracing approaches and improves computation speed even further with the same strict mathematical description.

METHODOLOGY
Our approach is based on the methodology described in Rofallski and Luhmann (2022). This approach shifts the ray tracing optimization problem for flat port housings from image space -requiring additional complex iterative procedures -to P p' u p' * A n-1 N 1 t 1 (X d 1 ,r d 1 ) µ 1 µ 2 µ n-1 µ n N n-1 A 2 A 1 p' m c t n-1 (X d n-1 ,r d n-1 ) Figure 1. Forward ray tracing through n media for a single point. From the measured point p ′ m , distortion correction is applied to gain p ′ u . The outgoing ray ⃗ Ai (red) is then refracted at n − 1 interfaces ti, considering refractive indices µi and µi+1 and the normal vectors at their intersection points ⃗ Ni. An unrefracted ray (green) with the projected point p ′ * is depicted for reference. object space to overcome these limitations. Our approach extends these findings by including hemispherical dome ports, assuming two interfaces of spheres. The approach is generic to cover other cases with, e.g. multiple interfaces, simultaneous camera calibration or object-invariant interfaces.
For alignment of the camera objective with the dome port, the generic optical system representation has to be considered. This includes the description with entrance and exit pupil through which all rays on the object side (entrance) and image side (exit) pass, assuming a straight ray path before and after the pupils. In photogrammetry, the pinhole model is widely used and represents a simplification of the former. In this case, no lens in front of the aperture is present and hence perspective center and entrance pupil coincide (Luhmann et al., 2020). In general, the entrance pupil and the dome center have to be aligned to obtain a refraction-free image which is the basis for the physical alignment. In a mathematical representation with a single viewpoint however (which is the case for the following), deviations from the perfect alignment are expressed with respect to the mathematical perspective center.

Ray Tracing
As shown in Figure 1, a measured image coordinate p ′ m is first undistorted at given approximate or calibrated values, according to its (approximate) interior orientation parameters in air (IOP). The model for the IOP can be any model that allows for undistortion of image coordinates. We used the common model by Brown (1971). The procedure is performed iteratively using Gauß-Newton optimization for correction of radial-symmetric and tangential-asymmetric (decentering) distortion, affinity and shear and principal point shift (Luhmann et al., 2020). The undistorted point p ′ u is then transformed to a vector from the given sensor position passing through the perspective center by employing parameters of interior and exterior orientation. This creates an arbitrarily scaled directional 3D vector ⃗ A1.
With a given starting point P0 i = (X0 i , Y0 i , Z0 i ) T , the image ray can be represented as a vector in the global coordinate system. For the first medium, the approximate position of the perspective center O ′ is used whereas the intersection points of the ray with the next interface is used for all subsequent media: (2) For n media, the dome port is parametrized as n − 1 concentric spheres (i.e. centers ⃗ Figure 1). The values are expressed bundle-invariant in the local camera coordinate system w.r.t. the perspective center. The fixed radii can usually be obtained from the manufacturer's data sheet. The searched intersection point ⃗ P0 i+1 on a sphere, given by its center ⃗ X d i and radius r d i in the local coordinate system, fulfills the following relation: The local dome offset coordinates ⃗ X d i are transformed to the global coordinate system by a 6DOF Helmert-Transformation. All subsequent Equations are related to the global coordinate system. Inserting Equation 2 into Equation 3 yields the following quadratic equation which has to be solved for the ray's scaling parameter si (Hanrahan, 1989): Obviously, this is a quadratic equation of the form as 2 i + bsi + c = 0 and can be solved analytically: As the origin of the ray lies inside the sphere, two possible solutions can be obtained from the positive definite discriminant. Since we are only focused on the points in front of the camera, a positive value for si has to be obtained for determination of the intersection point. Since all vectors from the sphere center to a point on the sphere intersect the surface perpendicularly, the normal direction at each intersection point ⃗ P0 i+1 is given by the vector from the sphere center to the point: With these parameters, we refract the 3D vector at the first spherical interface which results in the directional vector that travels through the next medium (Glassner, 1989): This relation is valid for all interface types that allow for calculating a normal vector and an intersection point. The algorithm from Equation 2 to 7 is performed recursively until all involved media are passed, e.g. twice for the common bundle-invariant case, air-glass-water. The last vector with a point on the outer interface ⃗ P0 n represents the outgoing ray An (Figures 1 and 2). This algorithm is performed for all exterior orientations (EOP) k and observed object points (OP) j.

Error function in object space
As stated in Rofallski and Luhmann (2022), the error function for optimization (e.g. for bundle adjustment) can be formulated in object space ( Figure 2). Hence, the orthogonal distance between the outgoing ray ⃗ An and the approximate value for the corresponding object point is obtained according to Mulsow and Maas (2014) and minimized to obtain a 3D residual vector ⃗ ∆xyz k,j , as follows: Pj j th point in object space ⃗ P0 n (k, j) Starting point of the outgoing ray for point j and exterior orientation k ⃗ An Direction vector of outgoing ray Contrary to usual bundle adjustment error functions where the reprojection error is minimized in image space, this has a handful of advantages, as pointed out by Rofallski and Luhmann (2022). Most importantly, after computationally cheaper undistortion, no iterative model coordinate has to be obtained in image space which is a major improvement over other ray tracing procedures in terms of computation speed and analytical evaluation, such as in Mulsow and Maas (2014).
The model was integrated into our high-performance bundle adjustment program, based on the C++ library Ceres-solver (Agarwal et al., 2022). In further contrast to Rofallski and Luhmann (2022), we integrated the model with the use of analytic derivatives for all parts, except for the iterative undistortion of the measured image coordinate. Strict use of the chain rule, as is the case with "automatic differentiation" in Ceres, obtains the strict derivatives with respect to all unknown parameters, resulting in a significant speed improvement over numerically determined (i.e. finite differences) derivatives. However, the undistortion part of the computation either fills in seamlessly with the analytical derivatives, if distortion parameters are optimized or they can be fully omitted if Visualization of the error function in object space. The orthogonal line-point distance ∆xyz k,j is calculated as a 3D vector that is subject to minimization by means of least squares. distortion parameters are not estimated in the optimization process. This results in an even higher computation speed gain.
After obtaining the resulting Jacobian J and the residual vector ⃗ dl, the normal equation system can be formulated: J Jacobian containing partial derivatives P Weight matrix ⃗ dx Reduced vector of unknowns ⃗ dl Reduced vector of observations For datum definition, we use Helmert constraints for translation and rotation and distance constraints for scale definition and hence a free network adjustment. The constraints are formalized in the matrix of constraints B and the vector of discrepancies ⃗ w. These are bordered on the normal equation matrix and absolute term vector (Luhmann et al., 2020).
For computational stability, we employ an implementation of the well-known Levenberg-Marquardt algorithm with multiplicative extension, according to Marquardt (1963). A pseudo code for the entire method is shown in Algorithm 1.

ERROR-FREE SIMULATION
To motivate for the evaluated data sets, we first show the theoretical effects of dome decentering in a simulation with error-free image coordinates. These were obtained by employing the ray tracing model by Mulsow and Maas (2014) and explicitly projecting the refracted image coordinate to the sensor. This procedure was performed iteratively, based on a small close-range data set with a GSD of about 30 µm and an average image scale of 1 : 5 (section 5.1). The major interest here is the resulting residual pattern when neglecting refraction.
We have simulated three data sets with a dome center shift of each 5 mm for each direction and an additional data Algorithm 1 Pseudo code for the optimization algorithm to optimize the squared sum of weighted residuals SSR set with all three coordinate components shifted by 5 mm. Both modeled spheres coincided, i.e. ⃗ X d 1 = ⃗ X d 2 . Apart from the camera constant set to c = −11 mm, all parameters were simulated with zero value. Subsequently, the data sets were adjusted using the Brown model with radial-symmetric and tangential-asymmetric (decentering) distortion, as well as affinity and shear to show the capability of absorbing the refraction parameters. Figure 3 shows the median residuals in a predefined grid over all images. Quite prominently, lateral shifts in X and Y cause systematic effects in the order of 1/3 px. In contrast, the longitudinal deviation in Z does not cause any major effects although small systematic patterns persist in the edges. The highest residual distortions arise from a shift in all coordinate directions up to 1 px. Table 1 shows the accuracy and precision details, comparing to the simulated reference coordinates. Relative accuracy is determined as the relation of RM SXY Z with the maximum object extent of 106.4 mm. Principally, a   perfectly axis-symmetrical behavior would be anticipated for shifts in δX and δY and point-symmetrical for δZ which is only approximately the case. The residual asymmetries most likely arise from an irregular point and camera station network from the original data set which was the base for this simulation.
Results in object space confirm the systematic residuals in image space and show a considerable accuracy loss which cannot be fully absorbed, especially for the laterally shifted data sets. The object space residuals (exemplary shown in Figure 4 for δX = δY = δZ = 5 mm) after transformation show a systematic wave-like radial pattern with outside vectors pointing in and inside vectors pointing out. Hence theoretically, errors arise from neglecting refraction in the camera model, Figure 4. Object space residuals compared to reference data set for δX = δY = δZ = 5 mm especially if lateral shifts are present. This however is on the margin that most likely only affects measurements of highest quality.
It has to be considered that the simulated image coordinates are free of any errors and might scale with measurement uncertainty in image space or the dome offset in the respective coordinate direction. However, increased noise in the data might as well "smear" the shown patterns and conceal the systematic effects that occur with smaller margin. To investigate on this behavior, real calibration data sets are evaluated in the next section.

EVALUATION ON REAL DATA SETS
Two data sets with different dimensions were evaluated to show possible limitations and applications for the model with real data. First, a close-range data set with a small flat calibration artifact (60 mm × 60 mm) was evaluated. Secondly, we focused on a medium-size calibration with a more three-dimensional calibration artifact of considerably larger size (1500 mm × 1000 mm × 200 mm). The camera systems and exemplary images are provided in Figure 5.
In addition to our approach, we also implemented the approach, given by Mulsow and Maas (2014) for hemispherical dome ports. Since the underlying physics are equal to our representation and only the optimization procedure is changed from image space to object space, it is assumed that results deviate only in small margins. However, the changed optimization procedure constrains the derivatives to be numerical and hence should require significantly more time per iteration than our model. The denomination for our approach will be "Ray tracing with minimization in object space" (RTO) whereas the standard by Mulsow and Maas (2014) will be referred to as "Standard ray tracing" (RT).
For all evaluations, we set the refractive index of air to µair = 1.00028 and for glass to µ glass = 1.49. Additionally, the dome radii were fixed, as stated from the data sheets of the respective dome ports and presented in Table 2. Apart from that, all parameters (i.e. exterior orientations [EOP], interior orientation [IOP], object points [OP], position of coinciding dome centers and refractive index of water) were part of adjustment if not stated otherwise for a given data set.
For precision verification, we considered internal statistics, such as standard deviation of unit weight σ0 and respective standard deviations for some parameters, obtained from the cofactor matrix. Accuracy was evaluated differently between the data sets. In data set 1, a set of reference coordinates was available. Resulting coordinates of each data sets were transformed by Helmert transformation without scale (6DOF)  Figure 5. Camera systems of the two data sets (a, b) and an exemplary image from the bundles (c, d). Note: Due to poor focus, the elevated points visible in (c) for the close-range data set were removed, resulting in the described 2D artifact.
to the reference coordinates, and RMS values of residuals were analyzed. For the second data set, 12 accurately determined reference lengths were available that were not used as constraints for adjustment. These were evaluated towards the length measurement error (LME) between actual and reference lengths, according to German guideline VDI/VDE 2634-1. Stated RMS and MAX values refer to the LME. The average uncertainty of the reference lengths was 2.8 µm in air. Not all of the stated reference lengths by Menna et al. (2018b) were available, hence LME errors are not fully comparable to the mentioned contribution.

Close-range calibration (data set 1)
This data set was acquired in a small water tank, filled with clear fresh water and originally obtained for calibrating an ultra-close-range underwater weld-inspection system by Kahmen and Luhmann (2022). A single monochromatic camera with 10 mm lens observed a flat calibration artifact that was marked with ring-coded targets. It was integrated in a low-cost acrylic housing with a 3" dome port which was not customized for the used camera and its placement was determined iteratively by the procedure in Kahmen and Luhmann (2022). Hence, offsets might be increased, compared to a commercial and customized dome housing that is available for DSLR cameras. The GSD was about 30 µm.
Placed about 100 mm below the water surface, the artifact was imaged at an average distance of 50 mm following a standard self-calibration protocol with convergent images and roll diversity which resulted in 12 images. Due to the limited depth of field in such close range, camera tilting could only coarsely be performed, as targets quickly run out of focus. Additionally, illumination was very challenging which led to underexposing targets in image corners. Three well-spread length constraints were introduced for scaling, as stated best practice by Luhmann et al. (2020). All object points observed with less than four rays and exterior orientations with less than six corresponding points in object space were excluded from the bundle. A dry calibration in air was obtained before the dome port was mounted onto the acrylic tube and the calibration underwater was performed. The IOP of the dry dataset were provided as fixed parameters for further analyses.

Medium-distance calibration (data set 2)
The presented data set resulted from a camera calibration, carried out at Baia del Rogiolo in Tuscany, Italy. For reference, the data was part of earlier contributions (Menna et al. 2017a;Menna et al. 2018b;Menna et al. 2018a).
A 1500 × 1000 mm 2 aluminum Dibond composite board with 160 circular coded targets was imaged with a Nikon D750 DSLR camera and 24 mm lens in a commercial dome port housing. This housing was customized for the specific camera-lens setup and should hence provide a well-centered configuration. The dome radius was more than twice the size of the radius from the first data set. Some elevated points at about 20 cm above the plane were present for additional depth variation in the artifact. The relevant parameters are summarized in Table 2. The test field was placed about five meters below water surface and imaged at an average distance of 1.2 m following a standard self-calibration protocol with convergent images and roll diversity which resulted in 25 images. The Ground Sample Distance (GSD) was 0.3 mm. Again, all coordinates with less than four image rays and images with less than six observed points were eliminated from the bundle beforehand. Figure 7 provides an overview of the imaging geometry.

Results
We analyzed the data sets with the following parameters: First, as a reference, we obtained the Brown model calibration with all camera and distortion parameters taking part in the optimization (BRN). Secondly, we optimized for the dome offsets while keeping the interior orientation constant either at values, provided by in-air calibration (RT/RTO Air) or those resulting from the aforementioned Brown calibration with the same data points as a two-step procedure (RT/RTO UW). Since no dry calibration was available for the second data set, the air-analyses were only performed for the first close-range data set. It is worth noting that between the in-air calibration (without dome) and the underwater calibration, the camera remained stable. Additionally, we also simultaneously optimized for both the dome offsets and the full camera parameters in the Brown model (RT/RTO + IOP). Results are depicted in Table 3. The analyses were performed on a new notebook with Intel i9, 10th Gen., 2.40 GHz, 64Gb RAM. No observations were weighed according to any apriori standard deviations to solely focus on the performance of the models.
For the close-range data set, an overall high precision and accuracy can be stated from σ0 and RM SXY Z in all analyses. The Brown approach had a σ0 value of about 1/12 px. The own implemented RTO algorithm produces similar results to the RT algorithm, apart from the values for the RTO + IOP estimation. These differ strongly between the two ray tracing algorithms. In general, estimating the dome offsets improves the accuracy, depending on the chosen configuration by 13 -34 %. The best performance was obtained from employing the two-step calibration by first solving for the Brown parameters and then fusing these as fixed parameters into the RT/RTO model. Calibration with simultaneous determination of IOP and the dome offsets were generally solvable in both models but produced highly unreliable results which differed significantly from all other parameter estimations. We observed that estimated parameters were, according to the parameters' standard deviations, significant but are highly doubtful concerning their physical meaning. Dome offsets were estimated significantly, except for the estimation with predetermined Brown parameters. In this case, offsets were rather small in lateral direction and at least one order of magnitude lower than with the in-air calibration for the longitudinal Z-direction. The processing time of the Brown model and the RTO model was equally fast at about 0.02 s without estimation of IOP and about factor ten slower when simultaneously estimating interior parameters, as well. Compared to standard ray tracing, this is a computation time decrease of about factor 80 and 40, respectively.
For the medium-size data set, high accuracy and precision persisted. The value of σ0 in image space was about 1/4 px.
Relative accuracy was about three times higher than with the close-range data set. However, no accuracy or precision gain was obtainable by introducing the dome offset parameters by RT and RTO. Dome offset parameters again were consistent among the ray tracing algorithms for the sole estimation of the offset despite not being estimated significantly. Adding the IOP to the problem causes unstable results and only rough accordance between the two ray tracing algorithms, as either the absolute value seems quite high (RT) or the standard deviation is very high (RTO). The mean LME is negative throughout all analyses at about 0.03 mm. Speed of the RTO algorithm was again equal to the Brown model when fixing interior orientation parameters. If these are optimized for as well, the time per iteration again decreases by factor ten in the RTO model. RT was again comparably slow and computation time relations from the first data set are reflected here, again.

Discussion
The two presented data sets are quite diverse in their dimensions, accuracies and known parameters. Firstly, we were able to determine offsets in all analyses but with different physical meanings and statistical reliability. Regarding physical interpretability, an in-air calibrated camera with added dome center should be calibrated to obtain meaningful parameters. Hence, dome offset parameters from RT/RTO (Air) are most likely closest to the actual parameters. This strict representation improved the accuracy of data set 1 but so did other analyses, as well. Unfortunately, we were unable to obtain a dry calibration for the second data set, disabling an independent verification. Next, calibrating the Brown parameters first with additional ray tracing is not strictly physical, for the calibrated distortion parameters absorb most of the resulting distortion. Residual offsets only add an additional radial correction term with a different principal point and small margin to account for any residual error patterns. This might as well cause overparametrization, especially in case of data set 2 where accuracy is decreased by introducing ray tracing parameters.
Results from the first and other (unpresented) data sets indicate an accuracy improvement from this procedure, despite its questionable physical meaning. Simultaneous estimation of IOP and dome offsets led to degenerated systems of equations. Major correlations of the dome offset with the principal point and camera constant arose and were close to |ri| = 1 (δX with x ′ p , δY with y ′ p and δZ with c). We tried excluding these parameters and adjusting for distortion parameters only but could not obtain any results of higher accuracy from this setting, either. This remained, even though all major correlations were eliminated with that step.
The ambiguous accuracies between data sets and accuracy loss for the ray tracing approaches in data set 2 may result from many factors which cannot be fully decided within this paper. A possible reason may be the degrading image quality towards the edges. As mentioned by Menna et al. (2018b), effects such as chromatic aberration, astigmatism and coma might cause an image quality loss which is significantly higher than refractive effects. Despite being a considerably bigger issue for flat ports, imperfectly centered dome ports suffer from these effects, as well. Also, the field curvature effect (i.e. a spherical focus field, enhanced by the lens character of the dome port) creates blurry areas in image corners (Menna et al., 2017b) which should cause the main image quality loss. Figure 8 shows a snippet of data set 2 from the image center and the top left corner. A strong decrease in sharpness can be observed in the top left image (b) whereas the center image contains sharp targets. Despite resulting from a customized housing and a high-quality DSLR camera, the sharpness decrease and refraction effects might overlay each other. It can thus be assumed that the centering is quite well and offsets are too small to be determined.
Additionally, the determinability might be affected by the given dome radii. With the first data set, rather small radii were present whereas the second data set had twice as large dome radii. Due to a higher curvature in a smaller radius, offsets might increase effects in the image and hence their determinability in the presented models. Also, the relation between the acquisition distance and the dome offset might have an effect.
The mean LME in the medium-size data set is negative throughout all analyses (except for the rather unreliable RTO + IOP), pointing towards a possible temperature-induced scaling effect. This might be caused by the water temperature affecting the calibrated lengths. However, temperature correction was performed beforehand by applying the given linear expansion coefficient of aluminum Dibond panels α = 0.024 mm/m · K to the lengths. However, temperature gradients can still occur and affect the assumed length deviations, probably causing the shown offset for the mean LME.
In the close-range data set, a comparably low relative accuracy is presented. This should be caused by the very large image scale which pushes towards the limits of photogrammetry. This concerns imaging quality, camera and object stability, heterogeneous illumination and even the targets' printing quality where very small errors cause a high relative error magnitude due to small object dimensions.
Computation time was significantly reduced with the RTO model. The algorithm can be described fully analytical if no interior orientation parameters are optimized. Hence, the method is as fast as the standard Brown model with the added value of obtaining a strictly physical model to the underwater analysis. Adding IOP to the problem, a numerical module is added which solves for the undistortion and its derivatives numerically and adds the result to the analytical values. Hence, a considerable speed decrease can be observed which still provides a speed improvement on the order of factor 40 -80 over the standard ray tracing. This speed improvement is even considerably higher than the stated factors by Rofallski and Luhmann (2022) who found a factor of roughly ten between their numerical implementation of the algorithm for flat ports and its counterpart by Mulsow and Maas (2014). Thus, analytical derivatives decrease computation time by factor eight in this case. We refrained from optimizing the refractive index as it is highly correlated and no useful results came from the inclusion of the parameter.

CONCLUSION
A versatile and efficient optimization model for dome ports, capable of providing statistical metrics in object space, was presented in this contribution. We have shown the theoretical relevance with simulated data. Furthermore, we presented two real data sets with different outcomes regarding accuracy improvement and parameter determinability. Computation time decreased significantly while results were similar to a standard ray tracing algorithm from literature. Thus, a significant improvement towards larger data sets and real-time applications is made for underwater analysis without loss of accuracy.
The presented computation time decrease can also be transferred to the flat port model by integrating analytical derivatives into the given model, as well. It is assumed that computation time would be comparable to the Brown model implementations, hence enabling real-time applications for both flat and dome ports with this improvement. Especially applications with ROV and other unmanned systems, relying on automated visual navigation (i.e. Simultaneous Localization and Mapping) might benefit from this model.
Large Structure from Motion configurations which include several thousands or millions of observations are currently impractical to solve with the standard ray tracing procedures.
Here, the proposed methodology might improve significantly towards reasonable computing times. We will focus on this in future work. Of special concern might be a trivial solution to the error function which hasn't raised any issues for us, so far. When coinciding all points and perspective centers in one single point, all rays pass exactly through all object points and residual errors become infinitely small. Only distance constraints or ground control points would cause any discrepancy in the model. Hence, good approximate values and probably more than one distance constraint or several ground control points are necessary to avoid this trivial solution, especially when many points are present. As stated before, this has not yet been an issue with any data set but should be kept in mind when designing the bundle geometry.
Fundamentally, the concentric sphere assumption, as well as fixed radii, can be eliminated if the bundle configuration allows for meaningful results. However, it has been shown that parameter determination is already rather complicated with small dome offsets. Hence, an even higher degree of freedom which would result from estimating additional radii or sphere centers may be unrewarding for future work. Especially changes in the radii have shown no observable accuracy influence and are very likely to not produce any meaningful results, if adjusted simultaneously. Additionally, the values can be gained from manufacturers' data sheets.
The accuracy gain in the presented data is rather small and probably not of practical concern. However, both our own presented simulations with error-free image coordinates and other simulations, e.g. by Kunz and Singh (2008), indicate a possible accuracy gain if the offset between entrance pupil (or perspective center in our model) and the dome center is high enough. A dependence on the dome radii might be possible and thus an explicit model be the more important for small radii. Simulation series on shifted dome offsets, different dome radii, measurement accuracy in image space and other functional constraints should be carried out in the future and can form a sharper picture on the matter.
The problematic determinability of the small offset parameters can complicate a closed control loop for automated camera centering inside dome port housings. Constraining one or two coordinate directions of the offset might be sufficient with conventional manufacturing accuracies. With only one remaining free direction of offset, results might stabilize and visual servoing be possible with the presented method.
All in all, we have presented a novel methodology with high potential towards large data sets which are currently hard to evaluate with the state-of-the-art ray tracing models. Dome offsets can be estimated from the evaluated data sets, further investigations towards the limits and advantages of such methodology have to be performed. However, the main advantage of the major computation time decrease persists and should benefit the multimedia photogrammetry community.