<figure><img src="image_1.png" /></figure>

HAL
open science

Refractive Pose Refinement

Xiao Hu, François Lauze, Kim Steenstrup Pedersen

► To cite this version:

Xiao Hu, François Lauze, Kim Steenstrup Pedersen. Refractive Pose Refinement. International Journal of Computer Vision, 2023, 131, pp.1448-1476. 10.1007/s11263-023-01763-4. hal-05526670

HAL Id: hal-05526670

https://hal.science/hal-05526670v1

Submitted on 25 Feb 2026

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# Refractive Pose Refinement

## Generalising the geometric relation between camera and refractive interface

Xiao Hu · François Lauze · Kim Steenstrup Pedersen

Received: date / Accepted: date

**Abstract** In this paper, we investigate absolute and relative pose estimation under refraction, which are essential problems for refractive structure from motion. To cope with refraction effects, we first formulate geometric constraints for establishing iterative algorithms to optimize absolute and relative pose. By classifying two scenarios according to the geometric relationship between the camera and refractive interface, we derive the corresponding solutions to solve the optimization problems efficiently. In the scenario where the geometry between the camera and refractive interface is fixed (e.g., underwater imaging), we also show that the refractive epipolar constraint for relative pose can be established as a summation of the classical essential matrix and two correction terms caused by refraction by using the virtual camera transformation. Thanks to its succinct form, the resulting refractive epipolar constraint can be efficiently optimized. We evaluate our proposed algorithms on synthetic data showing superior accuracy and computational efficiency compared to state-of-the-art (SOTA) methods. We further demonstrate the application of the proposed algorithms in refractive structure from motion on real data. Our datasets

(Hu et al., 2022b) and code (Hu et al., 2022a) are publicly available.

**Keywords** Refraction · Pose estimation · SfM · Reconstruction

## 1 Introduction

Camera pose estimation is one of the fundamental components for many important applications such as Structure from Motion (SfM), Visual Odometry (VO), and Multiview Stereo (MVS). During the past few years, we have witnessed tremendous successes achieved by SfM (Schonberger and Frahm, 2016) in homogeneous spaces (uniform media). However, less attention has been paid to the extension of SfM in non-uniform media, until recently. Light gets refracted when it passes through the interface of different transparent media, which makes it bend and change trajectories. This effect causes classical methods to malfunction (Treibitz et al., 2011). Examples for when this occurs include underwater imaging (Jordt et al., 2016), imaging through a water surface (Xiong and Heidrich, 2021a), or imaging objects encased in a transparent medium (Cassidy et al., 2020; Hu et al., 2021), as illustrated in Fig. 1. We also include illustrations of some use cases from the efforts of digitizing natural history collections (Hedrick et al., 2020), such as insect fossils included inside Baltic amber (Sadowski et al., 2021) and specimens in glass jar and alcohol. As refraction invalidates camera models (e.g., Fig. 2) and introduces systematic model errors, SfM needs to be adapted to it. In some early works, the effect is modeled as an error to perspective imaging and approximately corrected by focal length adjustment and radial lens distortion (Kang et al., 2012a; Lavest et al., 2000). Unfortunately, this solution only

<footnote>Xiao Hu
DTU Space, Technical University of Denmark & Department of Computer Science (DIKU), University of Copenhagen, Denmark
E-mail: hux062303@gmail.com</footnote>

<footnote>François Lauze
Department of Computer Science (DIKU), University of Copenhagen, Denmark
E-mail: francois@di.ku.dk</footnote>

<footnote>Kim Steenstrup Pedersen
Department of Computer Science (DIKU) & Natural History Museum of Denmark (NHMD), University of Copenhagen, Denmark
E-mail: kimstp@di.ku.dk</footnote>
<figure><img src="image_3.png" /><figcaption>**Fig. 1** Top row: (Left) Bathtub scene - 1 of 4 images used in the camera pose estimation and reconstruction in Fig. 28. (Right) A stag beetle encased in crystal resin - 1 of 30 images used for SfM reconstruction in Fig. 30. Bottom row: Two examples of applications of scenario 2 within digitization of natural history museum specimens. (Left) Insect fossil inclusion in Baltic amber submerged in syrup from the amber collection at Natural History Museum of Denmark (NHMD, 2022a). (Right) Specimens in glass jars and alcohol from the herpetology collection at Natural History Museum of Denmark (NHMD, 2022b).</figcaption></figure>

works in limited cases (Treibitz et al., 2011) due to the fact that the refractive distortion not only depends on the scene depth, but also the relative pose between the camera and the refractive interface. Recent work on the Refractive SfM (RSfM) problem tend to model the refractive effects explicitly and then adjust SfM pipelines accordingly, such as the absolute pose estimation approaches by Haner and Åström (2015) and Jordt and Koch (2013), the relative pose estimation solutions by Chadebecq et al. (2017, 2019) and Jordt and Koch (2013), and the adapted bundle adjustment solution by Jordt and Koch (2013). Applications of RSfM for underwater 3D reconstruction and endoscopy have been reported by Jordt et al. (2016); Chadebecq et al. (2017) and Chadebecq et al. (2019). Most of the previous works (Kang et al., 2012a; Lavest et al., 2000; Haner and Åström, 2015; Jordt and Koch, 2013; Chadebecq et al., 2017; Jordt et al., 2016; Chadebecq et al., 2019) are focused on the case of **planar interfaces** that is also studied in this work. Nonetheless, our proposed methods can be extended to other cases as long as the refractive interface is known.

Revisiting classical SfM pipelines results in that several components need to be adapted to refraction: 1) geometric verification in feature matching; 2) absolute and relative pose estimation; 3) triangulation; 4) bundle adjustment. As has been shown by Chari and Sturm (2009), classical essential, fundamental, and homogra-

<figure><img src="image_4.png" /><figcaption>**Fig. 2** The pinhole camera model and the corresponding variations under refraction.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Fig. 3** The considered scenarios: Scenario 1, the refractive interface is fixed in the camera coordinate frame and it moves together with the camera in the world coordinate frame, while in scenario 2, the refractive interface is fixed in the world coordinate frame.</figcaption></figure>

phy matrices do not hold under refraction, which means that, theoretically, geometric verification based on those matrices cannot be applied. We experimentally found that they can still be used, but at the expense of removing a certain number of correct matches. Triangulation needs to be computed using rays after refraction. The difficulty of applying bundle adjustment using forward projection lies in solving a very large amount of expensive quartic reprojection equations, which has been alleviated by Jordt and Koch (2013). Therefore, feature matching verification, triangulation and bundle adjustment could be considered as solvable to some extent. By contrast, absolute and relative pose estimation remains a challenging problem. Classical Perspective-N-Point (PnP) algorithms for absolute pose estimation, and the eight-point (Hartley and Zisserman, 2004) and five-point algorithms (Nistér, 2004) for relative pose estimation cannot work due to the violation of the pinhole camera model. As shown by Agrawal et al. (2012), the imaging system under refraction becomes an axial camera that is a special case of the generalized camera (GC) model (Grossberg and Nayar, 2001; Sturm, 2005). Unfortunately, previously proposed pose estimation methods for the GC model have been shown to be quite sensitive to noise, and the camera pose cannot be accurately estimated even in the low noise case (Chadebecq et al., 2019; Jordt, 2013).

In this paper, we focus on absolute and relative pose estimation of a camera under refraction. Two scenarios can be distinguished. They are shown in Fig. 3 and are as follows:
1. Scenario 1: The refractive interface is fixed in the camera coordinate frame, e.g., underwater imaging system, see for instance Chadebecq et al. (2017, 2019); Jordt et al. (2016); Jordt and Koch (2013); Oleari et al. (2015); Rizzini et al. (2015).

2. Scenario 2: The refractive interface is fixed in the world coordinate frame, e.g. an object embedded in a transparent medium (resin, amber *etc*. ) or underwater objects imaged above the surface, see for instance Hu et al. (2021).

Most of the previous works (Chadebecq et al., 2017, 2019; Jordt et al., 2016; Jordt and Koch, 2013) target scenario 1, except Ichimaru et al. (2019) that target scenario 2 instead. Our work considers pose estimation in both scenarios. Our primary contributions are two iterative minimisation algorithms for absolute and relative pose estimation under refraction, which works in both scenarios. Their formulations allow for easy derivations of the corresponding Jacobian matrices. Moreover, for relative pose estimation, we also establish a compact epipolar constraint under refraction using the virtual camera transformation. Specifically, in scenario 1, the epipolar constraint can be represented using a refractive essential matrix that is a summation of the classical essential matrix and two correction terms due to refraction. Evaluations on synthetic data are used to assess the performance of the proposed algorithms against state-of-the-art methods for noise robustness, accuracy, and computational efficiency. Results show the proposed method outperforms state-of-the-art methods. Furthermore, we show their applications in RSfM on real data. Our methods are more stable and efficient than previous methods targeting scenario 1 (Chadebecq et al., 2017, 2019; Jordt et al., 2016; Jordt and Koch, 2013). Compared with Ichimaru et al. (2019) who targets only a very restricted subset of scenario 2 as the parameters of the refractive interface in each camera view are supposed to all be the same, equal to the parameters of the refractive interface of the reference view, our work does not rely on such an assumption. As an extension of our previous work (Hu et al., 2021), where we primarily addressed the pose estimation problem in scenario 2, this paper not only extends our pose estimation methods to scenario 1, but also presents several novel findings. Overall, the contributions and advances over our previous work include:

- We extend our previous work (Hu et al., 2021) to accommodate scenario 1, which, therefore, improve the versatility of the proposed solution so that it can be applied in more user cases such as underwater imaging system and fluid-immersed endoscopy (Chadebecq et al., 2019).

- An efficient, accurate and robust to noise absolute camera pose estimation method.

- A simpler (than Chadebecq et al. (2019)) formulation of the epipolar constraint under refraction. In scenario 1, we show that this formulation can be represented with a refractive essential matrix that is composed of a classical essential matrix and two refraction caused correction terms.

- Specially, in scenario 1, we derive initialization solutions for absolute and relative pose estimation, which can estimate camera pose correctly when the noise level is small.

- We evaluate our methods on more real datasets, both on third-party datasets and our collected real datasets. Moreover, we will make our own datasets (Hu et al., 2022b) for refractive structure from motion and an open-source implementation of our solutions (Hu et al., 2022a) publicly accessible.

Though we derive our algorithms by assuming a single planar refractive interface, the proposed methods are also accommodated in the case of multi-layer refraction or non-flat refraction as long as the refractive interfaces are known. Overall, we believe the proposed solution contributes to the refractive pose estimation problem and the further development of RSfM applications in e.g. oceanic engineering, digitization of natural history collections, and digitization of cultural heritage and art.

The paper is structured as follows: Section 2 discusses the related work and Section 3 presents notations and preliminaries used in this paper. The proposed algorithms are further detailed in Section 4. Then, experimental results are presented and analyzed in Section 5. Finally, conclusion and future work are summarized in Section 6.

## 2 Related work

In this section, we briefly review related work on absolute and relative pose estimation as well as RSfM.

**Absolute Pose estimation:** Absolute pose estimation for a single perspective camera is well studied with a number of PnP solutions that achieve impressive accuracy, e.g. minimal solutions by Kneip et al. (2011); Gao et al. (2003) and non minimal solutions by Lepetit et al. (2009); Hesch and Roumeliotis (2011); Ferraz et al. (2014). In the presence of refraction, the geometry of a perspective camera no longer holds, which makes PnP solutions fail to estimate camera pose correctly. A few works have been proposed to tackle this problem. E.g., an eight-point algorithm to calibrate the refractive interface and estimate the camera pose was introduced
in Agrawal et al. (2012). However, when the refractive interface is known, this prior knowledge cannot be used by their method to remove degrees of freedom (DoF). Using accelerometers to estimate the camera's vertical direction, a two-point algorithm to estimate translation heading angle was proposed by Chang and Chen (2011), but limited to a horizontal refractive plane. A five-point algorithm using co-planarity constraints (Agrawal et al., 2012) is presented by Haner and Åström (2015) and demonstrated good performance. Nonetheless, our experiments showed that it is susceptible to noise. An alternating optimization that refines pose and depth sequentially is proposed by Jordt and Koch (2013), while our experiments show this solution is very sensitive to initialisation. As has been shown by Agrawal et al. (2012), the camera under refraction in scenario 1 works as an axial camera, which is a special case of the generalized camera model (Grossberg and Naray, 2001). Therefore, theoretically, algorithms developed for generalized cameras or multi-camera systems can possibly be applied for absolute pose estimation in scenario 1. In the lens of pose estimation for generalized cameras, a few papers investigate the absolute pose estimation problem for generalized cameras or multi-camera systems; Nistér and Stewénius (2007); Kneip et al. (2013); Lee et al. (2016); Miraldo et al. (2018) and Ventura et al. (2014) separately proposed a minimal solver for recovering the pose of a multi-camera system. Besides minimal solvers, non-minimal solvers (Kneip et al., 2014, 2013; Sweeney et al., 2014; Fragoso et al., 2020) have been proposed by solving polynomials using Gröbner basis solvers (Kukelova et al., 2008). However, they cannot be applied in scenario 2. Despite the existence of a few pose estimation methods for generalized cameras, as we will show in Section 5, they are very sensitive to noise.

**Relative Pose Estimation:** Chari and Sturm (2009) studied the two-view refractive geometry and derived corresponding refractive fundamental (rF) and homography matrices. Although this work is of theoretical interest, no practical application has been given. A seven-point linear solution for relative pose estimation with known vertical direction from an accelerometer was proposed in Chang and Chen (2011). When the camera orientation is known, Kang et al. (2012b) solves the relative translation optimally under the $L_\infty$-norm, and otherwise use an evolutionary algorithm for hybrid optimization. In the case that the orientation is not known, the authors proposed to use an evolutionary algorithm for hybrid optimization, however this is very time consuming. Unlike previous methods that require extra equipment or prior information, Jordt and Koch (2013)

introduces an alternating optimization using the geometric constraints proposed by Agrawal et al. (2012). Their solution has been integrated into an RSfM for deep-sea 3D reconstruction (Jordt et al., 2016). Similar to their absolute alternating optimization, this method also relies on very good initialisation. Chadebecq et al. (2017) shows that the two-view relationship under refraction can be established using the generalized epipolar constraint (GEC) (Pless, 2003). Although the GEC problem has minimal solutions (Stewenius et al., 2005; Ventura et al., 2015) and non-minimal linear ones (Li et al., 2008; Mouragnon et al., 2009; Kneip and Li, 2014), they are particularly sensitive to noise (Chadebecq et al., 2017, 2019). Thus, Chadebecq et al. (2019) further proposes a novel formulation of the rF constraint, of size 21 × 12. Their experiments show that relative pose estimation based on the proposed rF constraint outperforms other state-of-the-art techniques. However, feature-dependency and high dimensionality make their method time-consuming.

**Refractive Structure from Motion:** RSfM has been studied in recent years, where the majority of approaches adapt standard SfM routines by considering the effects of refraction as distortions and then compensating for them through rectification. Adapting radial distortion and focal length accordingly can compensate for refraction to a certain extent (Lavest et al., 2000). However, it only works within a limited range and cannot completely exclude systematic biases, as has been reported by Kang et al. (2012b) and (Jordt and Koch, 2012). Additionally, previous work by Jordt and Koch (2013) have demonstrated that the errors become non-negligible when the refractive interface and the image plane are not strictly fronto-parallel. This has sparked the interest of researchers to develop RSfM approaches based on the ray-based model, such as the RSfM method with known vertical direction (Chang and Chen, 2011; Kang et al., 2012b). A technique that associates each 3D point to its perspective 2D feature through a virtual camera was proposed by Jordt and Koch (2013), which facilitates applying the non-linear refractive bundle adjustment without the need to compute complex forward projections. Recent works by Chadebecq et al. (2017, 2019) suggest further consideration of the rF constraint in the non-linear refinement, which has been shown to give superior performance. However, all these methods are targeting scenario 1 and are not applicable in scenario 2. Ichimaru et al. (2019) proposed a unified RSfM framework that claims to work with moving camera and moving refractive interface. However, their "static interface" targets only a very restricted subset of scenario 2, as the parameters of the refractive interface in each
camera view are supposed to all be the same, equal to the parameters of the reference view. In other words, this restricts the camera motions, which is reflected in their experimental setup where the camera moves only along a line/guiding rail. By contrast, our work does not rely on such an assumption.

## 3 Notations and Preliminaries

Scalars, vectors, and matrices are denoted by lowercase, bold lowercase, and bold uppercase letters, respectively. The identity matrix is denoted **I**. The *i*-coordinate of an object **p** is denoted by **p**<sub>*i*</sub>, and the list of *i*, *j* and *k* coordinates by **p**<sub>*i*,j,k*</sub>. The cross product is represented by $\times$. $SO(3)$ is the group of rotations. Its tangent space at **I** is the Lie algebra $so(3)$ of skew-symmetric matrices and is isomorphic to $R^3$, we denote it by $\vee$ and its converse by $\wedge$, following the convention by Ma et al. (2004). Our formulations make systematic use of Plücker coordinates (Sturm, 2005). As the following derivations involve multiple coordinate frames and transformations, we introduce, for clarity, new notations:

- By default, W represents the world coordinate frame, C represents the camera coordinate frame, and V represents a virtual camera coordinate frame.

- In general, $^bR_a$ denotes a rotation matrix from coordinate frame *a* to coordinate frame *b*. $^bC_a$ denotes the origin of coordinate frame *a* in coordinate frame *b*. By default, **R** denotes a rotation matrix from the world frame to the camera frame (i.e. $^C R_W$) and **c** denotes the location of the camera center with respect to the world frame (i.e. $^W C_C$). $R_V$ denotes a rotation matrix from a virtual camera frame to its camera frame (i.e. $^C R_V$) and $c_V$ denotes the location of the virtual camera center with respect to the corresponding camera frame (i.e. $^C c_V$).

- $\mathbf{a}_1^i$ represent a variable (e.g. a point or a ray) indexed by *i* and defined in coordinate frame 1.

### 3.1 Background on Refraction

Snell's law, $\mu_i \sin \theta_i = \mu_j \sin \theta_j$, describes the relationship between the angles of incidence $\theta_i$ and refraction $\theta_j$ from the normal **n** and the refractive indices $\mu_i$ and $\mu_j$ of the media. Setting $\lambda = \mu_i/\mu_j$, Snell's law in vector form is (Fig. 4):

$$ \mathbf{r} := S(\lambda, \mathbf{n}, \mathbf{i}) = \lambda \mathbf{i} + \mathbf{n} \sqrt{1 - \lambda^2 [1 - (\mathbf{n} \cdot \mathbf{i})^2]} - \lambda(\mathbf{n} \cdot \mathbf{i})\mathbf{n}, \quad (3.1) $$

where **i**, **r** are the direction of the incident and transmitted ray, respectively. We have actually switched the

<figure><img src="image_6.png" /><figcaption>**Fig. 4** Illustration of Snell's law in vector form (left) and the idea of virtual camera coordinate (Right). The virtual camera transformation seeks a virtual camera that satisfies the classical perspective projection model.</figcaption></figure>

incidence and refraction direction but the "Snell" function above satisfies the symmetry relation

$$ \mathbf{r} = S(\lambda, \mathbf{n}, \mathbf{i}) \iff \mathbf{i} = S(\lambda^{-1}, \mathbf{n}, \mathbf{r}). \quad (3.2) $$

Furthermore, along with the origin **c** of the incident ray and the last interface parameter **d** (any point **p** on the plane satisfies $\mathbf{n}^\top\mathbf{p} + \mathbf{d} = 0$), the point of refraction **q** can be computed as

$$ \mathbf{q} = \left[ \left( \mathbf{c} + \frac{|\mathbf{n}^\top \mathbf{c} + \mathbf{d}|}{\mathbf{i} \cdot \mathbf{n}} \mathbf{i} \right)^\top \quad 1 \right]^\top. \quad (3.3) $$

These relations can be cascaded to model flat multilayer refraction, see Agrawal et al. (2012). In this work, we focus on the single-layer case though our proposed method can also be used for multiple-layer refraction as long as the parameters of the refractive interfaces are known.

### 3.2 Virtual Cameras Under Refraction

Virtual camera coordinates were first introduced by Telem and Filin (2010) and later extended by Jordt and Koch (2012). One can find a virtual camera frame such that the correspondence between each 3D point and its 2D feature is given by a perspective projection, see (Jordt and Koch, 2012) and Fig. 4. The virtual camera center $c_V$ is the intersection point of the refraction ray **r** and the line along the interface normal **n** through the real camera center **c**

$$ c_V = c + \alpha_2 n = q + \alpha_1 r. \quad (3.4) $$

The rotation matrix $R_V$ from the virtual camera coordinate frame to the camera coordinate frame is found by aligning the interface's normal to the optical axis (i.e., $e_3 = [0, 0, 1]^\top$):

$$ R_V = \exp(\hat{\mathbf{a}}), \quad \mathbf{a} = \theta \mathbf{l}, \quad \theta = \mathbf{e}_3 \cdot \mathbf{n}, \quad \mathbf{l} = \frac{\mathbf{e}_3 \times \mathbf{n}}{\|\mathbf{e}_3 \times \mathbf{n}\|}. \quad (3.5) $$
$\mathbf{R}_V$ is feature-independent; it only relies on the normal of the interface and the optical axis. $\mathbf{c}_V$ is however feature-dependent, it also depends on $\mathbf{r}$. The virtual focal length $f_V$ is the distance $d$ between the camera center and the refractive plane (Jordt and Koch, 2013). The relation between a 3D point $\mathbf{p}_V$ in virtual camera coordinates and its projection $\mathbf{u}_V$ becomes

$$ \mathbf{u}_V = f_V \left[ \frac{\mathbf{p}_{V,1}}{\mathbf{p}_{V,3}} \frac{\mathbf{p}_{V,2}}{\mathbf{p}_{V,3}} \right]^{\top}. \quad (3.6) $$

The main advantage of using virtual camera coordinates is that the forward projection can be done without solving a quartic equations for finding refractive points. This becomes significant in projection-demanding tasks, e.g. bundle adjustment (Jordt and Koch, 2013), as massively solving quartic equations is very time-consuming.

## 4 Methods

In the rest of this work, we assume a single planar refractive interface with known planar parameters. For comprehension, we first present the constraints used for absolute and relative pose refinement geometrically. Then, we detail the exact formulations for the two scenarios.

### 4.1 Absolute Pose Refinement

The geometric constraint shown in Fig. 6 is used to establish an optimization objective for absolute pose refinement. Given a 3D point $\mathbf{p}$ and its corresponding image point $\mathbf{u}$, if the parameters of the refractive interface are known, the line $l$ can be found. Clearly, $\mathbf{p}$ should lie on $l$, this will serve as the geometric constraint for absolute pose refinement. We now express the constraint using the Plücker coordinates for the line $l$. If the refractive ray $\mathbf{r}$ and the refractive point $\mathbf{q}$ are known, the line $l$ parallel to the ray direction $\mathbf{r}$ and passing the refractive point $\mathbf{q}$ is given by its Plücker coordinates as:

$$ l = \begin{bmatrix} \mathbf{q}_4 \mathbf{r}_{1,2,3} - \mathbf{r}_4 \mathbf{q}_{1,2,3} \\ \mathbf{q}_{1,2,3} \times \mathbf{r}_{1,2,3} \end{bmatrix} = \begin{bmatrix} \mathbf{q}_4 \mathbf{r} \\ \mathbf{q}_{1,2,3} \times \mathbf{r}_{1,2,3} \end{bmatrix} =: \begin{bmatrix} l_1 \\ l_2 \end{bmatrix} \quad (4.1) $$

since we can assume $\mathbf{q}_4 \neq 0$ ($\mathbf{q}$ is in the Euclidean space) and $r_4 = 0$ (it's a direction). The 3D point $\mathbf{p}$ lies on $l$ if and only if

$$ l_1 \times \mathbf{p} + l_2 = \mathbf{0}. \quad (4.2) $$

We will use a least-squares formulation of this constraint for the two scenarios, but in scenario adapted frames:

<figure><img src="image_7.png" /><figcaption>Fig. 5 Illustration of the geometric constraint for absolute pose refinement in two scenarios. We highlight the pose-dependent variables using red color.</figcaption></figure>

- In scenario 1, as shown in Fig. 5, the refractive interface is known in the camera frame. Therefore, by (3.2), the direction $\mathbf{r}$ and the point of refraction $\mathbf{q}$ could be computed using feature's coordinates and the refractive interface in the camera frame. Further, $l$ can be calculated using the direction $\mathbf{r}$ and the point of refraction $\mathbf{q}$. Apparently, $l$ in the camera frame does not depend on the camera pose. We thus choose to represent (4.2) in the camera frame. Under this representation, it is $\mathbf{p}$ that depends on the camera pose ($\mathbf{R}, \mathbf{c}$), $\mathbf{p} = \mathbf{R}(\mathbf{p}_W - \mathbf{c})$.

- In scenario 2, the refractive interface is known in the world frame (i.e., $\mathbf{n}$ and $d$ are known in the world frame). We therefore represent (4.2) in the world frame. The computation of $l$ in the world frame is a bit involved. We first have the interface normal in camera frame $\mathbf{n}_C$ computed as $\mathbf{n}_C = \mathbf{R}^{\top}\mathbf{n}$. Then, for a given incidence direction $\mathbf{i}_C$ in the camera frame, we get the refracted ray $\mathbf{r}_C$ as $\mathbf{r}_C = S(\lambda, \mathbf{R}^{\top}\mathbf{n}, \mathbf{i}_C)$. Additionally, the equation of the refractive interface is $\mathbf{n}_C^{\top}\mathbf{x}_C + d_C = 0$ with $d_C = d + \mathbf{n}^{\top}\mathbf{c}$. By (3.3), the refractive point $\mathbf{q}$ becomes:

$$ \mathbf{q}_C = \left[ \frac{|d + \mathbf{n}^{\top}\mathbf{c}|}{\mathbf{i}_C^{\top}\mathbf{R}^{\top}\mathbf{n}} \mathbf{i}_C^{\top} \ 1 \right]^{\top}. \quad (4.3) $$

As $\mathbf{r}_C$ and $\mathbf{q}_C$ have dependencies on the camera pose, so does $l$. Compared with scenario 1, the involvement of camera pose in (4.2) is more complicated due to the non-linearity of Snell's law.

Regardless of the scenario-adapted formulations, with (4.2), this provides a generic least-square objective function of ($\mathbf{R}, \mathbf{c}$) for both scenarios:

$$ \arg\min_{\mathbf{R}, \mathbf{c}} \sum_{i=1}^{N} \| \mathbf{l}_1^i \times \mathbf{p}^i + \mathbf{l}_2^i \|^2 = \sum_{i=1}^{N} \| \mathcal{L}^i(\mathbf{R}, \mathbf{c}) \|^2, \quad (4.4) $$

with $\mathcal{L}^i(\mathbf{R}, \mathbf{c}) = \mathbf{l}_1^i \times \mathbf{p} + \mathbf{l}_2^i$, its dependency in ($\mathbf{R}, \mathbf{c}$) being described above, while $i$ being the index of an image point and $N$ being the total number of image points.
Refractive interface

<figure><img src="image_8.png" /><figcaption>**Fig. 6** Illustration of the geometric constraint for absolute pose refinement. The point **p** lies on the line **l**, which is used to formulate a cost function for absolute pose refinement.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Fig. 7** The epipolar constraint under refraction using the corresponding virtual cameras. For an arbitrary feature point in view 1 and its corresponding match in view 2, two virtual cameras are established. The refraction rays (r₁ and r₂) and the ray connecting two virtual cameras' centers (r₃) are co-planar, which establishes the epipolar constraint under refraction.</figcaption></figure>

Note that as (4.2) is given up to a multiplicative constant; such a constant need to be chosen in (4.4), which amounts to fixing **q**₄. There will be a natural choice for each scenario. Unlike Ichimaru et al. (2019) who minimize the reprojection error or the distance between 3D points, the objective in (4.4) is a sum of square distances from points to a given line.

## 4.2 Relative Pose Refinement

In the following, we first introduce an epipolar constraint and then the optimization objective using the epipolar constraint for relative pose refinement.

### 4.2.1 Epipolar Constraint from a New Perspective

Unlike previous works (Chari and Sturm, 2009; Chadebecq et al., 2017, 2019; Jordt and Koch, 2013; Jordt et al., 2016) that derive an epipolar constraint using the Klein quadric, we establish the epipolar constraint by

extending the classical epipolar constraint using the virtual camera transformation. As shown in Fig. 7, given an arbitrary feature point in view 1 and its matched feature point in view 2, two corresponding virtual cameras are obtained. Furthermore, the rays after refraction (r₁ and r₂) and the ray connecting the two virtual cameras centers (r₃) are co-planar. If r₁, r₂, and r₃ are represented in a common coordinate frame, they satisfy the epipolar relation

$$ \mathbf{r}_2^\top \mathbf{r}_3^\top \mathbf{r}_1 = 0. \qquad (4.5) $$

### 4.2.2 Optimization Objective

Inspired by Chadebecq et al. (2019), we formulate a joint optimisation problem on the relative pose and tri-angled points by minimising the reprojection errors and the epipolar costs

$$ \mathbf{R}, \mathbf{c}, \mathbf{p}^i = \arg \min_{\mathbf{R}, \mathbf{c}, \mathbf{p}^i} \sum_{i=1}^{n} \mathcal{E}^i(\mathbf{R}, \mathbf{c}, \mathbf{u}^i, \mathbf{u}^{i'})^2 \\ + \sum_{i=1}^{n} \|\mathcal{R}^i(\mathbf{I}, \mathbf{0}, \mathbf{p}^i, \mathbf{u}^i)\|_2^2 + \|\mathcal{R}^i(\mathbf{R}, \mathbf{c}, \mathbf{p}^i, \mathbf{u}^{i'})\|_2^2, \qquad (4.6) $$

with **u**^i, **u'^i** being i-th feature point in view 1 and its match in view 2, respectively. $\mathcal{E}(\mathbf{R}, \mathbf{c}, \mathbf{u}^i, \mathbf{u}^{i'})$ represents the epipolar cost defined in the left hand side of (4.5) and the terms $\mathcal{R}(\cdot, \cdot, \mathbf{p}^i, \mathbf{u}^i)$ represent the reprojection error of the i-th point in the corresponding view.

For relative pose estimation, we assume that the first camera coordinate frame coincides with the world frame, i.e. its pose is $\mathbf{R}_W = \mathbf{I}$ and $\mathbf{W}\mathbf{c}_C = \mathbf{0}$. Unlike Chadebecq et al. (2019), which computes reprojection errors using forward projection by solving quartic equations, we compute reprojection errors via the virtual camera coordinates. Since we also formulate the epipolar costs using virtual camera coordinates, the inclusion of epipolar costs will not bring a heavy computational burden. By contrast, in Chadebecq et al. (2019), more computations are required to lift the coordinates. Compared with Ichimaru et al. (2019) who merely minimize reprojection errors, (4.6) jointly minimizes epipolar costs and reprojection errors, offering two improvements:

- The epipolar cost only relies on the relative pose (no dependency on the scene) so that it helps to optimize the relative pose regardless of the 3D points.

- As shown by Hadfield et al. (2018), using a combination of two cost surfaces may decrease the possibility of getting trapped at a local minimum, since it requires two cost functions to reach the same local minimum point at the same time.
## 4.3 Scenario 1

### 4.3.1 Absolute Pose Estimation

In this scenario, our optimisation problem, expressed in camera frame, is written

$$ \arg \min_{\mathbf{R},\mathbf{c}} \sum_{i=1}^{N} \|\hat{\mathbf{l}}_{1}^{i} \times \mathbf{R} (\mathbf{p}_{\mathrm{W}}^{i} - \mathbf{c}) + \mathbf{l}_{2}^{i}\|^{2} = \sum_{i=1}^{N} \|\mathcal{L}^{i}(\mathbf{R}, \mathbf{c})\|^{2}. \quad (4.7) $$

Here we set $\mathbf{q}_4 = 1$. This is a nonlinear least squares formulation on $\mathbb{S}^3(3) \times \mathbb{R}^3$. To solve it we use a manifold-adapted Levenberg-Marquardt algorithm (subsection 4.5). To do so, we need an expression for the Jacobian of $\mathcal{L}^i$ for each feature point $i$. A variation around $(\mathbf{R}, \mathbf{c})$ is written as $(\mathbf{R}e^{i\hat{\xi}}, \mathbf{c} + t\mathbf{v})$ with $\hat{\xi} \in \mathfrak{s}\mathfrak{o}(3)$. With it, the expression is

$$ J\mathcal{L}^i(\mathbf{R}, \mathbf{c})[\hat{\xi}, \mathbf{v}]^\top = -[\hat{\mathbf{l}}_1^i \mathbf{R} (\mathbf{p}_W^i - \mathbf{c})^\wedge, \hat{\mathbf{l}}_1^i \mathbf{R}][\hat{\xi}, \mathbf{v}]^\top, \quad (4.8) $$

i.e., $J\mathcal{L}^i$, in this system is represented by,

$$ J\mathcal{L}^i(\mathbf{R}, \mathbf{c}) = -[\hat{\mathbf{l}}_1^i \mathbf{R} (\mathbf{p}_W^i - \mathbf{c})^\wedge, \hat{\mathbf{l}}_1^i \mathbf{R}]. \quad (4.9) $$

Empirically, this optimization usually converges within a few iterations.

### 4.3.2 Initialization

As an iterative optimization, the refinement also needs an initialisation. It is known that, for non-convex optimization problems, the quality of the optimization result relies on the goodness of the initial value. A good initial value not only makes optimization converge faster but also drives the final result to a better local minimum. Therefore, it is crucial to develop an efficient and good initialization. In the following, we present a minor contribution, i.e. an efficient absolute pose estimation solver.

The proposed solver is inspired by famous EPnP by Lepetit et al. (2009) algorithm and can be treated as a variant of EPnP to handle the pose estimation under refraction. The EPnP (Lepetit et al., 2009) proposes using the barycentric coordinates to express a 3D point $\mathbf{p}$ as a weighted sum of four control points in world frame, i.e. $\mathbf{c}\mathbf{p}_W^j, j = 1, \dots, 4$:

$$ \mathbf{p} = \sum_{j=1}^{4} \alpha_j \mathbf{c} \mathbf{p}_W^j, \text{ where } \sum_{j=1}^{4} \alpha_j = 1. \quad (4.10) $$

<figure><img src="image_10.png" /><figcaption>**Fig. 8** Illustration of the Initialization solution for absolute pose estimation. Inspired by the EPnP algorithm that represents points with the help of four control points, four more control points are defined in the virtual camera frame so that points in the virtual camera frame can also be expressed using the barycentric coordinates.</figcaption></figure>

Next, by transforming $\mathbf{p}$ into the camera frame using camera pose $\mathbf{R}, \mathbf{c}$, we have the transformed point $\mathbf{p}_C$ as:

$$ \begin{aligned} \mathbf{p}_C &= \mathbf{R}^\top (\mathbf{p} - \mathbf{c}) = \sum_{j=1}^{4} \alpha_j \mathbf{R}^\top (\mathbf{c}\mathbf{p}_W^j - \mathbf{c}) \\ &= \sum_{j=1}^{4} \alpha_j \mathbf{c}\mathbf{p}_C^j, \text{ where } \mathbf{c}\mathbf{p}_C^j = \mathbf{R}^\top (\mathbf{c}\mathbf{p}_W^j - \mathbf{c}). \end{aligned} \quad (4.11) $$

By coordinate lifting, the EPnP algorithm cleverly encapsulated camera pose into the lifted coordinates of the four control points and finally formulate the absolute pose estimation problem as a linear system of the form $\mathbf{M}\mathbf{x} = \mathbf{0}$.

In scenario 1, to cope with refraction, we propose to further define four control points in the virtual camera frame $\mathbf{c}\mathbf{p}_V^j$ as:

$$ \mathbf{c}\mathbf{p}_V^j = \mathbf{R}_V^\top \mathbf{c}\mathbf{p}_C^j, \quad (4.12) $$

where $\mathbf{R}_V, \mathbf{c}_V$ be the transformation of the corresponding virtual camera coordinate system with respect to the camera coordinate system. Note that, in the transformation, we do not include $\mathbf{c}_V$ since we have shown in Section 3.2 that $\mathbf{c}_V$ is feature-dependent. Finally, as shown in Fig. 8, for any point $\mathbf{p}_V$ in the virtual camera frame, we have:

$$ \mathbf{p}_V = \mathbf{R}_V^\top (\mathbf{R}^\top (\mathbf{p} - \mathbf{c}) - \mathbf{c}_V) = \sum_{j=1}^{4} \alpha_j \mathbf{c}\mathbf{p}_V^j - \mathbf{R}_V^\top \mathbf{c}_V \quad (4.13) $$
Recall that by using the virtual camera coordinate (3.6), the refraction ray and $\mathbf{p}_V$ satisfy the perspective projection model:

$$z \begin{bmatrix} \mathbf{u}_V \\ 1 \end{bmatrix} = \mathbf{K}_V \mathbf{p}_V, \qquad (4.14)$$

where $z$ is the feature depth and $\mathbf{K}_V$ is the intrinsic matrix of the established virtual camera. By (4.13), the normalized coordinates (i.e. $[\tilde{\mathbf{u}}_V, 1]^\top = \mathbf{K}_V^{-1} [\mathbf{u}_V, 1]^\top$) and the Kronecker product, we can rewrite (4.14) as

$$\underbrace{\alpha \otimes \begin{bmatrix} 1 & 0 & -\tilde{\mathbf{u}}_{V_1} \\ 0 & 1 & -\tilde{\mathbf{u}}_{V_2} \end{bmatrix}}_{\mathbf{A}} \mathbf{x} = \underbrace{\begin{bmatrix} \mathbf{t}_{V_3} \tilde{\mathbf{u}}_{V_1} - \mathbf{t}_{V_1} \\ \mathbf{t}_{V_3} \tilde{\mathbf{u}}_{V_2} - \mathbf{t}_{V_2} \end{bmatrix}}_{\mathbf{b}}, \quad (4.15)$$

where $\mathbf{t}_V = -\mathbf{R}_V^\top \mathbf{c}_V$, $\mathbf{x} = [\mathbf{c}_V^{1\top}, \dots, \mathbf{c}_V^{4\top}]^\top$, and $\alpha = [\alpha_1, \alpha_2, \alpha_3, \alpha_4]^\top$. It is worthy to mention that, in scenario 1, as the parameters of the refractive interface is known in the camera frame, $\mathbb{R}_V$ and $\mathbf{v}$ are known and can be precomputed. Thus, $\mathbf{A}$ and $\mathbf{b}$ are known and decoupled with $\mathbf{x}$, which results in a linear system of $\mathbf{x}$ to be solved. By establishing (4.15) for each 3D-2D correspondence and stacking the results together (assume $N$ correspondences in total, $N \ge 6$), we will have a stacked linear equation:

$$M\mathbf{x} = \mathbf{h}, \qquad (4.16)$$

where $M = [\mathbf{A}^{1\top}, \dots, \mathbf{A}^{N\top}]^\top$ and $\mathbf{h} = [\mathbf{b}^{1\top}, \dots, \mathbf{b}^{N\top}]^\top$. If $M$ is of full rank, $\mathbf{x}$ can be solved by the pseudo-inverse. We denote the solution using the pseudo-inverse as $\mathbf{x}^* = M^\dagger \mathbf{h}$.

<figure><img src="image_11.png" /><figcaption>**Fig. 9** Illustration of the optimization process for (4.17). The weights are optimized by aligning the control points recovered from the pseudo-inverse solution and null vectors $\mathbf{x} = \mathbf{x}^* + \sum_{j=1}^{N'} \lambda_i \mathbf{v}_i$, and the control points in world frame.</figcaption></figure>

However, in practice, $M$ might be rank deficient, e.g. in the co-planar case, the rank will be 9 in the noise-free condition. In the case that $M$ is not full rank, the solution will be the weighted sum of $\mathbf{x}^*$ and a varying number of the null vectors in the right null space of $M$, i.e. $\mathbf{x} = \mathbf{x}^* + \sum_{j=1}^{N'} \lambda_i \mathbf{v}_i$ where $N'$ is the dimension of

the null space, $\mathbf{v}_i$ are the null vectors, and $\lambda_i$ are the unknown weights. Inspired by Ferraz et al. (2014), we estimate the optimal solution by solving the following optimization problem:

$$\mathbf{x} = \mathbf{V}\beta^*, \quad \text{where: } \beta^* = \arg\min_{\beta} \|\mathbf{V}\beta - \hat{\mathbf{x}}\|^2, \quad (4.17)$$

where $\mathbf{V}$ is a $12 \times (N' + 1)$ matrix composed column-by-column of $\mathbf{x}^*$ and $\mathbf{v}_i$, $i = 1, \dots, N'$, $\beta$ is a $(N' + 1)$-dimensional vector of weights, $\hat{\mathbf{x}} = [\hat{\mathbf{c}}_V^{1\top}, \dots, \hat{\mathbf{c}}_V^{4\top}]^\top$ is the transformed control points in the virtual camera frame using the recovered pose. More specifically, as illustrated in Fig. 10, we first solve for $\mathbf{x}$ which are the four control points in the virtual camera frame. Next, we recover $\mathbf{R}$ and $\mathbf{c}$ by solving a 3D-3D registration problem (Arun et al., 1987). We then get the transformed control points by $\hat{\mathbf{c}}_V^j = \mathbf{R}_V^\top \mathbf{R}_V^\top (\mathbf{c}_w^j - \mathbf{c})$. The optimal estimate will minimize the residuals between the sets of control points. It is worthy to note that, once $\mathbf{V}$ has been established, the complexity of this optimization is independent of the number of features considered. Therefore, this solution is of the complexity $O(N)$.

## Relationship with Generalized PnP

- Interesting readers may have noticed that by using the virtual camera transformation, we actually establish a multi-camera system in which each feature point will be associated with a virtual camera. Thus, theoretically, generalized PnP solvers for multi-camera system (e.g. UPnP (Kneip et al., 2014)) can also be applied to solve the refractive PnP problem. However, we argue that the refractive PnP problem is more challenging than the PnP problem for multi-camera system. As can be seen in Nistér and Stewénius (2007); Lee et al. (2016); Miraldo et al. (2018); Kneip et al. (2014, 2013); Ventura et al. (2014), all generalized PnP solvers assume the relative translations and orientations from each camera to the common body frame are known with high accuracy. In the case of refractive PnP problem, this is equivalent to assume that the $\mathbf{R}_V$ and $\mathbf{c}_V$ for each feature are known accurately, or in other words, features are known accurately. Apparently, this assumption is not possible in practice because noises are inevitable so that the noises on the features will directly affect the computations of $\mathbf{R}_V$ and $\mathbf{c}_V$.

- The proposed initialization method appears analogous to the gPnP proposed by Kneip et al. (2013) for multi-camera systems. However, the gPnP recovers the solution by solving a polynomial system using the Gröbner basis, which is different with our
solution. Moreover, experimental results show that our solution is more accurate and robust than gPnP. we also found that gPnP cannot give correct results in the co-planar case even in a noise-free condition. By contrast, our solution also works in the co-planar case, see section 5.1.1.

### 4.3.3 Relative Pose Estimation

Before giving exact solutions to (4.6), we first write explicitly the refractive epipolar constraint in scenario 1. As shown in Fig. 7, assuming the transformation from the virtual camera frame to the real camera frame are given as $R_{V_i}, c_{V_i}$ ($i$ denotes the $i^{th}$ view) and the transformation from the view 2 to the view 1 is given as $R, c$, the three rays $r_1, r_2, r_3$ in the first virtual camera frame are computed as follows:

$$ 
\begin{aligned}
\mathbf{r}_1 &= \mathbf{R}_{V_1}^{\top} \mathbf{u}, \quad \mathbf{r}_2 = \mathbf{R}_{V_1}^{\top} \mathbf{R} \mathbf{u}' \\
\mathbf{r}_3 &= \mathbf{R}_{V_1}^{\top} (\mathbf{R} \mathbf{c}_{V_2} + \mathbf{c} - \mathbf{c}_{V_1}).
\end{aligned}
\quad (4.18) $$

Then, the co-planar constraint can be written as

$$ \mathbf{u}'^{\top} \mathbf{R}^{\top} \mathbf{R}_{V_1} (\mathbf{R}_{V_1}^{\top} (\mathbf{R} \mathbf{c}_{V_2} + \mathbf{c} - \mathbf{c}_{V_1}))^{\top} \mathbf{R}_{V_1}^{\top} \mathbf{u} = 0. \quad (4.19) $$

Since $(\mathbf{R}\mathbf{a})^{\wedge} = \mathbf{R}(\mathbf{a})^{\wedge} \mathbf{R}^{\top}$, we could simplify the constraint (4.19) as

$$ \mathbf{u}'^{\top} \underbrace{\mathbf{R}^{\top} (\mathbf{R} \mathbf{c}_{V_2} + \mathbf{c} - \mathbf{c}_{V_1})^{\top}}_{\mathbf{r}\mathbf{E}} \mathbf{u} = 0. \quad (4.20) $$

Informally, we call $\mathbf{r}\mathbf{E}$ the refractive essential matrix

<figure><img src="image_12.png" /><figcaption>Fig. 10 Geometric illustration of the refractive epipolar constraint (4.20) in Scenario 1. Due to refraction, the virtual camera center shifts from the camera center by $c_{V_1}$ and $c_{V_2}$, respectively. By considering these shifts, the refractive epipolar constraint can be established using coplanarity.</figcaption></figure>

in the following.

- Recall that, in the case of no refraction, the classical essential matrix is $\mathbf{R}^{\top}(\mathbf{c})^{\wedge}$. By comparing the two essential matrices, we could find that the refractive essential matrix is composed of the classical essential matrix and two correction terms due to refraction. If there is no refraction, (4.20) will return back to the classical epipolar constraint since

$\mathbf{R}_{V_i} = \mathbf{I}$ and $c_{V_i} = \mathbf{0}$. Note that this not only holds in the single layer flat refraction case, but also in more general cases (e.g. multi-layer flat refraction or non-flat refractive interface) as long as the parameters of the refractive interfaces are known.

- Compared with the SOTA work by Chadebecq et al. (2019), our formulation is of lower dimension ($3 \times 3$ vs $21 \times 12$). Thanks to its simple form, we will be able to derive the analytical Jacobian matrices which are very useful to boost the iterative refinement. As we will show in Section 5.1.2, our solution is almost 100 times faster than their method, but still achieves similar performance.

- The derived refractive essential matrix shares a similar problem with Chadebecq et al. (2019); it is feature-dependent. Therefore, to evaluate the epipolar cost, we need to re-compute the refractive essential matrix for each feature pair. Luckily, it is easy to see that the computation can be done in parallel.

With (4.20), we could derive a quadratic approximation of the objective in the local tangent space. The optimization objective is the summation of three error terms, i.e. the epipolar cost and two reprojection errors. Thus, it is only necessary to separately derive the residual function and the corresponding Jacobian matrix for the three terms. The derivation is similar to the derivation for absolute pose refinement. Thus, for brevity, we omit the details, and only show the final result. The residual function and the corresponding Jacobian matrix for the epipolar cost and the reprojection error are given in (A.1) and (A.2), respectively. Note that here we have $\mathbf{x} = [\boldsymbol{\xi}^{\top}, \mathbf{v}^{\top}, \mathbf{w}^{i\top}]^{\top}$, $i = 1, \dots, N$ which denotes a variation around $(\mathbf{R}, \mathbf{c}, \mathbf{p}^{i})$.

### 4.3.4 Initialization

By using (4.20), we present another minor contribution, i.e. a relative pose initialisation solver. It can estimate relative pose correctly in low noise conditions. First, we expand the epipolar constraint as follows

$$ \mathbf{u}'^{\top} (\mathbf{c}_{v_2}^{\wedge} \mathbf{R}^{\top} + \mathbf{R}^{\top} \mathbf{c}^{\wedge} - \mathbf{R}^{\top} \mathbf{c}_{v_1}^{\wedge}) \mathbf{u} = 0, \quad (4.21) $$

By replacing $\mathbf{R}^{\top}$ with $\mathbf{R}'$ and $\mathbf{c} = -\mathbf{R}'^{\top} \mathbf{t}$, we further have:

$$ 
\begin{aligned}
& \mathbf{u}'^{\top} (\mathbf{c}_{v_2}^{\wedge} \mathbf{R}' - \underbrace{\mathbf{t}^{\wedge} \mathbf{R}'}_{\mathbf{E}} - \mathbf{R}' \mathbf{c}_{v_1}^{\wedge}) \mathbf{u} = 0 \Rightarrow \\
& \underbrace{\left[ \mathbf{u}'^{\top} \mathbf{c}_{v_2}^{\wedge} (\mathbf{u}^{\top} \otimes \mathbf{I}) - \mathbf{u}'^{\top} ((\mathbf{c}_{v_1}^{\wedge} \mathbf{u})^{\top} \otimes \mathbf{I}) - \mathbf{R}_2^{\top} (\mathbf{u}^{\top} \otimes \mathbf{I}) \right]}_{\mathbf{A}} \\
& \left[ \begin{matrix} \mathrm{vec}(\mathbf{R}) \\ \mathrm{vec}(\mathbf{E}) \end{matrix} \right] = 0,
\end{aligned}
$$
<figure><img src="image_13.png" /><figcaption>**Fig. 11** Configuration of scenario 2. The known variables are drawn in black color, while the unknown variables are depicted in red color. As **n** and **d** are only known in the world coordinate frame, the refractive point **q** and the refraction ray **r** have dependencies on the camera pose **R**, **c** and cannot be computed without knowing them. By aligning the normal of the refractive interface to the gravity direction, we define a new world coordinate frame, which can help us simplify the following derivations.</figcaption></figure>

$$ (4.22) $$

where $\text{vec}(\mathbf{A})$ represents the vectorization operation that stacks matrix elements column by column.

To solve (4.22), we suggest using the method proposed by Li et al. (2008). Assuming $\mathbf{A}_E$ and $\mathbf{A}_R$ are sub-matrices of $\mathbf{A}$ consisting of the first and last 9 columns, we could find the solution that satisfies $||\text{vec}(\mathbf{E})|| = 1$ by solving

$$ (\mathbf{A}_R\mathbf{A}_R^+ - \mathbf{I}) \mathbf{A}_E\text{vec}(\mathbf{E}) = 0 \quad (4.23) $$

where $\mathbf{A}_R^+$ is the pseudo-inverse of $\mathbf{A}_R$. After having $\mathbf{E}$, we can decompose it to get $\mathbf{R}$, $\mathbf{t}$ (up to a unknown scale) and finally $\mathbf{c}$. If the noise is negligible, it is possible to recover the scale by substituting $\mathbf{R}$ and solving a least square problem to get $\mathbf{t}$, as has been reported in Jordt and Koch (2013) and Chadebecq et al. (2017).

## 4.4 Scenario 2

In scenario 2, the parameters of the refractive interface is known in world coordinates. In standard SfM, the world coordinate frame is usually taken as the first camera coordinate frame. To simplify a bit our calculations, we use another world coordinate frame, adapted to the refractive interface: We choose the frame such that the interface goes through the origin and its normal is $\mathbf{n} = \mathbf{e}_3$, as the $\mathcal{F}_I$ illustrated in Fig. 11. Once the refinement is performed; results can be transferred back to the original world frame. The transformation from the new world frame to the original world frame can be computed by first aligning $\mathbf{n}_W$ to $\mathbf{e}_3$ using $^W\mathbf{R}_I$

(computed by (3.5)) and then shifting the refractive plane so that the origin lies on it:

$$ ^W\mathbf{c}_I = -d_W \mathbf{n}_W, \quad (4.24) $$

where $\mathbf{n}_W$ and $d_W$ are the plane parameters in the original world frame. In the following, the world frame is designated to the aforementioned new world frame for brevity.

### 4.4.1 Absolute Pose Estimation

Let $\mathbf{K}$ be the camera intrinsic parameter matrix. The back-projected ray of an arbitrary image point $\mathbf{u}$ (in homogeneous coordinate) in the world frame is given by $\mathbf{i}_W = \mathbf{R}^\top \mathbf{i}_C$, where $\mathbf{i}_C = \frac{\mathbf{K}^{-1} \mathbf{u}}{\|\mathbf{K}^{-1} \mathbf{u}\|}$. In the world frame, computing the refractive point and the refraction ray is easily done using (3.1) and (3.3) (we can assume, without loss of generality, that $^W\mathbf{c}_C(3) < 0$, i.e. the camera is below the refractive plane)

$$ \begin{aligned} \mathbf{q}_W &= \begin{bmatrix} ^W\mathbf{c}_C \mathbf{i}_{W,3} - ^W\mathbf{c}_C \mathbf{i}_W \\ \mathbf{i}_{W,3} \end{bmatrix} \\ \mathbf{r}_W &= \begin{bmatrix} \lambda \mathbf{i}_{W,1} \\ \lambda \mathbf{i}_{W,2} \\ \gamma \\ 0 \end{bmatrix}, \end{aligned} \quad (4.25) $$

where $\gamma = \sqrt{(1 - \lambda^2) + \lambda^2 \mathbf{i}_{W,3}^2}$. The Plücker coordinate $\mathbf{l}_W$ of the line passing by $\mathbf{q}_W$ and parallel to $\mathbf{r}_W$ is obtained by (4.1).

Finally, the residual function of the objective defined in (4.4) for absolute pose estimation is given as

$$ \mathcal{L}(\mathbf{R}, \mathbf{c}) = \mathbf{l}_{W_1} \times \mathbf{p}_W + \mathbf{l}_{W_2}. \quad (4.26) $$

Note that the camera pose $\mathbf{R}, \mathbf{c}$ are now implicitly embedded in $\mathbf{l}_W$. Its Jacobian matrix $J\mathcal{L}(\mathbf{R}, \mathbf{c}) = [J\mathcal{L}_R, J\mathcal{L}_c]$ is given in Appendix A.2.1.

**Initialization** In scenario 2, there is no efficient method for initialisation, neither can we derive an efficient solution as we do for scenario 1. Therefore, we use the 8-point algorithm proposed by Agrawal et al. (2012) as an initialization solution.

### 4.4.2 Relative Pose Estimation

For relative pose estimation, after obtaining the inhomogeneous coordinates of the refractive point and the refraction ray in the world frame (i.e. $\mathbf{r}_W^1$ and $\mathbf{r}_W^2$) from (4.25), we compute their virtual camera coordinates following Section 3.2. Since $\mathbf{n}_W$ is $\mathbf{e}_3$ in the world
coordinate frame, the resulting rotation matrix $^W\mathbf{R}_V$ is the identity matrix $\mathbf{I}$. To find the virtual camera location $^W\mathbf{c}_V$, we need first to solve the following equation

$$ \mathbf{q}_W + \alpha_1 \mathbf{r}_W = \mathbf{w} \mathbf{c}_C + \alpha_2 \mathbf{n}_W . \qquad (4.27) $$

This equation has a closed-form solution, with $\alpha_1 = \frac{w\mathbf{c}_{C,3}}{\lambda_{W,3}}$. Therefore, the virtual camera's location $^W\mathbf{c}_V$ is

$$ \mathbf{w}\mathbf{c}_V = \left[ \mathbf{w}\mathbf{c}_{C,1} \quad \mathbf{w}\mathbf{c}_{C,2} \quad \frac{w\mathbf{c}_{C,3}}{\lambda_{W,3}} \gamma \right]^{\top} . \qquad (4.28) $$

Then the ray $\mathbf{r}_W^i$ that connects the two virtual camera centers is $\mathbf{r}_W^3 = \mathbf{w}\mathbf{c}_{V_2} - \mathbf{w}\mathbf{c}_{V_1}$. With $\mathbf{r}_W^1, \mathbf{r}_W^2$, and $\mathbf{r}_W^3$, we can now write down the corresponding residual functions.

**Epipolar Cost:** The residual function for the epipolar cost term is given by:

$$ \mathcal{E} = \mathbf{r}_W^2 \top (\mathbf{w}\mathbf{c}_{V_1} - \mathbf{w}\mathbf{c}_{V_1})^{\top} \mathbf{r}_W^1, \qquad (4.29) $$

where $\mathbf{w}\mathbf{c}_{V_i}$ denotes the virtual camera location of the $i$th view and $\mathbf{r}_W^i$ represents the refraction ray computed from the same view.

**Reprojection Error:** The reprojection residual is:

$$ \mathcal{R} = \begin{bmatrix} f_V \frac{\mathbf{p}_{V,1}}{\mathbf{p}_{V,3}} - f_V \frac{\mathbf{r}_{V,1}}{\mathbf{r}_{V,3}} \\ f_V \frac{\mathbf{p}_{V,2}}{\mathbf{p}_{V,3}} - f_V \frac{\mathbf{r}_{V,2}}{\mathbf{r}_{V,3}} \end{bmatrix} . \qquad (4.30) $$

The Jacobian matrices $\mathcal{J}\mathcal{E}_x$ and $\mathcal{J}\mathcal{R}_x$ are provided in Appendix A.2.2.

**Initialization:** There is no efficient and robust method for estimating the relative pose under refraction. A solution can be to simply ignore the refraction and use the classical 5-point algorithm proposed by Nistér (2004) to get the relative pose. Afterward, 3D points are triangulated using the refraction rays.

## 4.5 Optimization on Manifolds

The optimisations problems (4.4) and (4.6) are unconstrained optimization problems on manifolds of the type $(\mathcal{S}\mathcal{O}(3) \times \mathbb{R}^3)^n$, $n \ge 1$. We use the Riemannian Levenberg-Marquardt algorithm of Absil et al. (2009) tailored to this situation. Using rectified coordinates, i.e. translating the $\mathcal{S}\mathcal{O}(3)$-Jacobian components to $\mathcal{S}\mathcal{O}(3)$, the tangent space of $\mathcal{S}\mathcal{O}(3)$ at $\mathbf{I}$, via the matrix exponential and logarithm (Ma et al., 2004), and projecting and translating the increment back on the manifold, see Fig. 12 for a graphical illustration.

<figure><img src="image_14.png" /><figcaption>**Fig. 12** Left: Mapping relationship among Lie group, Lie algebra, and the vector space. Right: An illustration of optimization on manifolds.</figcaption></figure>

## 5 Experiments

In this section, we present the experimental results for both scenario 1 and scenario 2. For each scenario, we first compare the proposed solutions to SOTA absolute and relative pose solutions using numerical simulations. We evaluate the performance in terms of accuracy, robustness to noise, and computational efficiency. Next, we report the experimental results obtained on real data, where we compare poses obtained using the proposed method to the ground truth poses. Furthermore, we demonstrate on real data, the application of the proposed methods to reconstruct objects under refraction. All experiments were run on a 2.8 GHz Intel Core 2 Duo computer.

Throughout this section, we define the rotation error and translation error between the ground truth $\mathbf{R}_{gt}$, $\mathbf{t}_{gt}$ and the estimation $\hat{\mathbf{R}}$, $\hat{\mathbf{t}}$ as

$$ \varepsilon_{\text{rot}} = \| \log(\mathbf{R}_{\text{gt}}^{\top} \hat{\mathbf{R}} - \mathbf{I}) \|_{\text{fro}}, \quad \varepsilon_{\text{tran}} = \| \mathbf{t}_{\text{gt}} - \hat{\mathbf{t}} \|_{\text{fl}}. \qquad (5.1) $$

### 5.1 Scenario 1: Numerical Experiments

First we consider results based on numerical simulations for Scenario 1.

#### 5.1.1 Absolute Pose Estimation

We have chosen the following SOTA algorithms for comparison with our proposed method: 1) gPnP (Kneip et al., 2013); 2) gP+s (Ventura et al., 2014); 3) 5pt (Haner and Åström, 2015); 4) AGW (Agrawal et al., 2012); 5) MMPPE (Miraldo et al., 2018). We name the proposed initialization and refinement as PL, PL+O respectively. Note that we intended to include other SOTA solvers like UPnP (Kneip et al., 2014). However, it is found that UPnP does not perform well in the case of refraction even without noise addition, see Fig. 35 in Appendix B.2. For gPnP, gP+s, and 5pt, we use their
open-source implementations provided by Kneip and Furgale (2014); gP+s (2014); Haner and Åström (2015). For MMPPE, we integrate a RANSAC (Fischler and Bolles, 1981) preprocessing procedure to their implementation (MMPPE, 2018) in order to select the best estimation, while for AGW, we re-implement their method to tackle the single layer refraction case. Except for gPnP that is implemented in C++, all the other methods are implemented in MATLAB.

**Non-Planar Case:** In the non-planar case, we conduct stochastic experiments by assuming a perspective camera with a focal length of 4800, a principal point of (960, 540), and no distortions. For each iteration, 100 uniformly distributed features are generated and perturbed by different levels of zero-mean Gaussian noise. We assume a flat refractive plane located at the origin with the normal being **n** = [0, 0, 1]ᵀ. To ensure the plane can be viewed in front of the camera, we draw the optical axis as **z** = [λ₁, λ₂, -1]ᵀ, λ₁, λ₂ ∈ [-0.5, 0.5] and then compute the camera orientation **R** so that **Rn** = **z**. We sample the camera location from a uniform distribution in the space that is 1-1.5 above the refractive plane and compute the rays after refraction using Snell's law. Finally, we sample world points from the rays after refraction, see Fig. 33 in Appendix B.1 as an example. The simulated refractive index is 1 (air) and 1.5 (glass) and we assume the camera is in the air.

In order to evaluate the noise resilience, we gradually increase the standard deviation of the Gaussian noise from 0 to 2 pixels. For each noise level, 100 simulations are executed. The resulting box plots show the rotation and translation error in Fig. 13, which also illustrate the mean runtime. As can be seen, our solution PL+O outperforms other SOTA solvers by showing the best noise resilience in all situations. Especially compared with the other iterative solver AGW, PL+O not only shows superior robustness, but also spends less computational time. It is worth noticing that AGW is intended to jointly calibrate the refractive interface (i.e. the normal and plane distance) and the camera pose. Its versatility results in considering unnecessary degrees of freedom if the parameters of the refractive interface are already known. Additionally, the proposed initialization PL also demonstrates similar performance to other competing methods.

**Planar Case:** In the planar case, we use the same setup for data generation, except we force all world points to be located on a common plane, see Fig. 34 for an example. Again, the noise level is varied from 0 to 2 pixels and 100 simulations are conducted per each noise level. The results are shown in Fig. 14. Note that we found

that gPnP, gP+s, and MMPPE do not work well in the planar case even without noise perturbation, see Fig. 36 in Appendix B.2. As indicated in Fig. 14, PL+O again maintains the best noise resilience and significantly outperforms another iterative alternative AGW. Regarding runtime, PL+O is much faster than AGW. The initialization solution PL shows good performance when the noise level is not more than 1 pixel, but degrades when the noise level continues to increase. One possible reason could be the dimension of the null space may grow over 4 due to the noise, e.g. Kneip et al. (2013) empirically found that the dimension of the null space may grow up to 6 in some cases.

### 5.1.2 Relative Pose Estimation

To generate a synthetic dataset for evaluating relative pose solvers, we consider the following setup. We assume two cameras with similar configurations (focal length and resolution) to the one used in the absolute pose case. We assume the interface is parallel to the image plane. 100 features are uniformly sampled in the first view and their corresponding 3D points are uniformly distributed along the refractive ray at a distance between 5 and 10 meters. Next, we project 3D points into the second view to get the feature points, see Fig. 35 for an example. Zero mean Gaussian noise is sampled and added to feature points in both views. We vary the standard deviation of the noise from 0 to 2 pixels to evaluate the robustness to noise. For each noise level, 50 simulations are executed.

We compare our method with four SOTA methods: The five-point algorithm (5pt) (Nistér and Stewénius, 2007) by ignoring the refraction, the seventeen point algorithm (17pt) for solving the GEC problem (Li et al., 2008), another two recent iterative refinement solvers using GEC (Chad-GEC) (Chadebecq et al., 2017) and the refractive fundamental matrix (Chad-rF) (Chadebecq et al., 2019). We intended to add other GEC solvers such as Kneip and Li (2014) and Stewenius et al. (2005), however, they do not perform well even in the noise-free condition, see Fig. 37 in Appendix B.3. We name our method Virtual Epipolar Constraint – VEC – in the following. Regarding implementation, for the 17pt and 5pt algorithms, we use the implementations provided by the OpenGV (Kneip and Furgale, 2014) library and the OpenCV (2022) library. For the Chad-GEC and Chad-rF algorithms, we had to re-implement due to the lack of available open-source code. All other solvers are implemented in MATLAB. For all the iterative solvers, we initialize the relative pose using the 5pt algorithm by ignoring the refraction. World points are then initialized by triangulation. Previous works by Chadebecq et al.
<figure><img src="image_15.png" /><figcaption>**Fig. 13** Comparison of absolute pose solvers for scenario 1 with respect to varying noise levels in the nonplanar case: the left box plot shows rotation error; the middle box plot shows the translation error; the right figure shows the computational time of all algorithms. Note that the mean runtime of AGW is beyond the scope of the y-axis when the noise level is above zero.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>**Fig. 14** Comparison of absolute pose solvers for scenario 1 with respect to varying noise levels in the planar case: the left box plot shows rotation error; the middle box plot shows the translation error; the right figure shows the computational time of all algorithms. Note that the mean runtime of AGW is beyond the scope of the y-axis.</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>**Fig. 15** Comparison of relative pose solvers for scenario 1 with respect to varying noise levels: the left box plot shows rotation error; the middle box plot shows the translation error; the right figure shows the computational time of all algorithms. Note that 17pt gives very accurate pose estimation results in the noise-free condition so that its box is very tiny.</figcaption></figure>
<figure><img src="image_18.png" /><figcaption>**Fig. 16** Comparison of 3D trajectory on the underwater ArUco dataset (scenario 1) by Zhang et al. (2021). The red dashed line is the trajectory estimated using the method by Zhang et al. (2021), and the blue dashed line refer to estimated trajectory using our absolute pose estimation method PL+O. The green curve is the ground truth trajectory as provided with the dataset.</figcaption></figure>

(2019) and Jordt and Koch (2013) have shown that the scale can not be reliably recovered, except in extremely low noise situations. Therefore, we estimate the scale by $s = \sum_{i}^{N-1} \hat{d}^i/d^i$, where $d$ (and $\hat{d}$) is the Euclidean distance between two (estimated) points.

Experimental results are shown in Fig. 15. We observe several points from Fig. 15: (1) The 17pt GEC solver shows excellent performance in the noise-free case. However, its performance degrades dramatically with the increase of the noise level. This confirms the statement by Chadebecq et al. (2017) and Jordt and Koch (2013) that GEC is very sensitive to noise. (2) Except in the noise-free case, 5pt shows better performance compared to 17pt in estimating the rotation. That is also why we choose to use the 5pt algorithm for initialization. (3) The three iterative refinement methods all demonstrate improvements in optimizing the relative pose. The proposed VEC outperforms the other two alternatives in terms of noise resilience and computational efficiency. Despite VEC shows superior performance and efficiency, we do observe 8 failure cases out of 450 simulations. The Hessian matrix of those failure cases is near singular, which makes the optimized results nonsense. This leaves us some future work to further improve the numerical stability of our method.

## 5.2 Scenario 1: Experiments on Real Data

Here we show results for real data representing Scenario 1.

<figure><img src="image_19.png" /><figcaption>**Fig. 17** Plots of localization errors.</figcaption></figure>

### 5.2.1 Real Data: Underwater ArUco Dataset

We use the dataset provided by Zhang et al. (2021) which includes an immersed ArUco (Garrido-Jurado et al., 2016) fiducial marker captured by a stereo camera encapsulated in a waterproof housing. The resolution of the collected images are 640 × 480 pixels and the configurations (normal, air-to-glass thickness, and glass thickness) of the waterproof housing was provided by the authors. Here we extend our absolute pose estimation approach to handle the double flat layer refraction, i.e. air-glass-water. The ArUco fiducial marker has size 16 centimeters and is used for pose estimation using the method proposed by Zhang et al. (2021) and our solution (only left camera is used). The results are illustrated in Fig. 16 where the red dashed line plots the trajectory estimated by Zhang et al. (2021), the blue dashed line plots the trajectory obtained by our solution, and the green solid line represents the ground truth. Furthermore, we show the localization errors in Fig. 17. It is from Fig. 16 and Fig. 17 that our solution shows significant superiority to their method, with much lower localization errors. The mean localization error of our solution is 1.59 cm, compared to 5.70 cm of their solution. The possible reason could be that they compute the pose by aligning 3D points triangulated from the stereo vision camera, which involves additional noise caused by stereo calibration parameters as well as noise existing in the right images (we only use the left camera images).

## 5.3 Scenario 1: SfM Application

To evaluate the contribution of our proposed method in refractive structure from motion (RSfM), we build an incremental SfM pipeline by replacing classical ab-
solute and relative pose estimation methods with the proposed algorithms. For pose estimation, we detected SIFT features (Lowe, 2004) from the segmented refractive images. Feature correspondences were then established with feature matching, where erroneous matches were removed by a geometric verification. We apply the refractive bundle adjustment from Jordt and Koch (2013) to refine camera poses and triangulated points. To get a relatively dense reconstruction result, we use the dense optical flow algorithm proposed by Kroeger et al. (2016) to get dense matches and triangulate them to get the ultimate point cloud.

### 5.3.1 Real Data: Underwater Checkerboard Dataset

This dataset was collected by manually moving an immersed GoPro Hero 9 Black camera with narrow lens setting (f/2.5) and 5184×3888 pixels resolution in front of a checkerboard and capturing static images at several places, as shown by examples in Fig. 18. The checkerboard is drawn on a piece of flat wood and held under water by hand. We calibrated the GoPro camera in air to obtain its intrinsics and distortion parameters. To get the interface parameters, we further did an underwater calibration using the method proposed by Agrawal et al. (2012). It is noted that it is unclear whether the housing of the GoPro camera should be modeled as a single- or double-layer refraction model. Therefore, we tried both cases, but found out the single-layer model always gave us lower reprojection errors. As a result, we finally decided to use the single-layer refraction model. The calibrated refractive interface parameters are [-0.102, 0.025, 0.994, -0.0098]. The final reconstruction results and the reconstruction result by AliceVision (2022) are shown in Fig. 19. Note that COLMAP (Schönberger and Frahm, 2016) cannot converge on this dataset. As can be seen from the side view, our simple incremental RSfM gives a more smooth result with a slightly more flat surface compared to a reconstruction not taking refraction into account. This is also consistent with previous findings reported by Chadebecq et al. (2019). Quantitatively, we fit a plane to each reconstructed point cloud and use the Root Mean Square Error (RMSE) of plane fitting as the metric for evaluation. As the point cloud is up to scale, we use the iterative closest point (ICP) algorithm to align their scales before applying plane fitting. The RMSE obtained using the proposed method is 0.0092, being much smaller than the RMSE obtained using AliceVision (being 0.0152), which also validates our observation.

### 5.3.2 Real Data: Underwater Dinosaur Dataset

The second dataset was collected similarly except that we replace the checkerboard with a clay dinosaur head, as shown in Fig. 20. The reconstruction results are shown in Fig. 21 based on 13 images. Again, we tried to use COLMAP (Schönberger and Frahm, 2016) on this dataset, but it failed to converge. As can be seen, AliceVision (2022) showed surprisingly very good results even though it completely ignores the refraction caused by the waterproof camera housing, better than the results obtained using our incremental RSfM. We found the possible reasons as follows:

- The implemented incremental RSfM pipeline cannot be directly compared with the performance of the pipeline used in AliceVision. We found that AliceVision can generate a large number of sparse points after structure from motion (approx. 2000 points) while our pipeline only gets approximately 200 sparse points. Moreover, in theory, we could include our linear relative pose estimation solutions into a robust estimator like RANSAC for rejecting outliers. But in practice, as the relative pose estimation solution that is based on the GEC is very sensitive to noises Jordt and Koch (2013); Chadebecq et al. (2017, 2019), we find out it is not very stable to use it for outlier removal. Therefore, this limits the number of features used for initialize the feature map that will be used for camera registration. By contrast, AliceVision can handle outliers very well.

- Our dense matching strategy (dense optical flow) and incremental strategy constrains us in using more views to adjust the reconstructed structure properly. As optical flow relies on the assumption of consistent intensity, it cannot produce dense matches of good quality on views that have severe intensity variations. The incremental scheme means if one image fails, the sequences after it will not be used. These disadvantages limit us to being able to use only a few views for dense reconstruction. More specifically, we used 13 images for structure from motion, i.e. estimating camera pose and a sparse point cloud. But since the dense optical flow results for image pair (5-6) is very bad, the remaining images were not used for dense reconstruction as wrong matches will result in many false points. By contrast, AliceVision uses multi-view stereo that can benefit from possibly all views, which helps it properly refine the 3D structure. This leaves us a second future improvement, i.e. integrate our solution into a mature SfM framework.

- As can be seen from the refractive interface parameters, the interface is almost orthogonal to the optical
<figure><img src="image_20.png" /><figcaption>**Fig. 18** Representative raw and segmented refractive images of the real underwater checkerboard dataset used for 3D reconstruction.</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>**Fig. 19** SfM reconstructions on the underwater Checkerboard dataset: (Left) SfM reconstruction based on refractive SfM using our method for Scenario 2. (Right) SfM reconstruction based on perspective camera using *AliceVision* (2022).</figcaption></figure>

axis and the distance from the camera center to the interface is only 9.8 mm. Considering that the ratio of the refractive index between the air and water is not very big, and the dinosaur head is viewed in the center of the image and the size of the dinosaur head is relatively small, the distortions caused by refraction are not severe enough so that classical methods can still give good estimates. This finding coincides with the previous works by (Jordt and Koch, 2013), (Chadebecq et al., 2017), and (4.6), where the discrepancy between the result of underwater RSfM (more specifically, single-layer underwater refractive case) and that of a classical SfM is limited.

This leaves us a third future work, i.e. to investigate the tolerance of classical methods to the distortions caused by refraction in certain circumstances. This will help us better understand the benefits of explicitly consider refraction, especially in the underwater scenario.

## 5.4 Scenario 2: Numerical Experiments

Next we turn our attention to Scenario 2 and again start with results based on numerical simulations.

### 5.4.1 Absolute Pose Estimation

The SOTA algorithms for comparison for absolute pose estimation include the 5pt algorithm (Haner and Aström, 2015), the AGW algorithm (Agrawal et al., 2012), and the classical EPnP algorithm (Lepetit et al., 2009). Regarding the implementation, we use the open-source code provided by Haner and Aström (2015) for the 5pt algorithm, while for the EPnP algorithm, we use the implementation provided by the OpenCV (2022) library. We re-implement the AGW algorithm for handling single layer refraction as the original implementation provided by Agrawal et al. (2012) was aiming to handle the double layer refraction case. The proposed method is named PO for short. We use the 8-pt algorithm (Agrawal et al., 2012) for initialization.
<figure><img src="image_22.png" /><figcaption>**Fig. 20** Representative raw and segmented refractive images of the underwater dinosaur dataset used for 3D reconstruction.</figcaption></figure>

<figure><img src="image_23.png" /><figcaption>**Fig. 21** SfM reconstructions on the underwater Dinosaur dataset: (Left) SfM reconstruction based on refractive SfM using our method for Scenario 2. (Right) SfM reconstruction based on perspective camera using AliceVision (2022).</figcaption></figure>

**Non-Planar Case** In the non-planar case, we create stochastic experiments by assuming a perspective camera with a focal length of 4800, a principal point of (960, 540), and no distortions. For each iteration, 100 uniformly distributed features are generated and perturbed by different levels of zero-mean Gaussian noise. We assume a flat refractive plane located at the origin with the normal being **n** = [0, 0.5, 1]ᵀ. To ensure the plane can be viewed in front of the camera, we place the camera above the refractive interface and randomly draw the optical axis as **z** = [λ₁, λ₂, -1]ᵀ, λ₁, λ₂ ∈ [-0.5, 0.5]. The simulated refractive index is 1 (air) and 1.5 (glass) and we assume the camera is in the air.

In order to evaluate the robustness to noise, we gradually increase the standard deviation of Gaussian noise from 0 to 2 pixels. For each noise level, 100 simulations are executed. In Fig. 22, the box plots show the rotation and translation error as well as the mean runtime. As can be seen, our solution PO outperforms other SOTA solvers by showing the best robustness to noise in all situations. When compared to the other iterative solver, AGW, PO not only shows superior robustness but also at a lower computational time cost.

**Planar Case** In the planar case, we use the same setup for data generation, except we force all world points to be on a common plane. Again, the noise level is varied from 0 to 2 pixels and 100 simulations are conducted per each noise level. The results are shown in Fig. 23, and PO again maintains the best robustness to noise and significantly outperforms other methods. Regarding runtime, PO is faster than AGW, but slower than the EPnP and 5pt algorithm.

### 5.4.2 Relative Pose Estimation

To generate a synthetic dataset for evaluating relative pose solvers, we first create a random refractive interface. Then, 100 features are randomly generated from a uniform distribution in the first view, and their corresponding 3D points are randomly uniformly distributed
<figure><img src="image_24.png" /><figcaption>**Fig. 22** Comparison of absolute pose solvers for scenario 2 with respect to varying noise levels in the nonplanar case: the left box plot shows rotation error; the middle box plot shows the translation error; the right figure shows the computational time of all algorithms. Note that, for 5pt and EPnP, their errors are out of scope.</figcaption></figure>

<figure><img src="image_25.png" /><figcaption>**Fig. 23** Comparison of absolute pose solvers for scenario 2 with respect to varying noise levels in the planar case: the left box plot shows rotation error; the middle box plot shows the translation error; the right figure shows the computational time of all algorithms. Note that, for 5pt and EPnP, their errors are out of scope.</figcaption></figure>

<figure><img src="image_26.png" /><figcaption>**Fig. 24** Comparison of relative pose solvers for scenario 2 with respect to varying noise levels: the left box plot shows rotation error; the middle box plot shows the translation error; the right figure shows the computational time of all algorithms. Note that, for 5pt, its translation errors are out of scope.</figcaption></figure>
<table><caption>**Table 1** Quantitative Result on Rubik's cube dataset: average reprojection error in pixels. Non refers to the non-refractive method.</caption><thead><tr><th>Solver</th><th>Non</th><th>5-pt</th><th>PO</th><th>AGW</th><th>EPnP</th><th>Ichimaru</th></tr></thead><tbody><tr><th>Error</th><td>54.6</td><td>7.6</td><td>5.9</td><td>266.9</td><td>407.2</td><td>611.7</td></tr></tbody></table>

on the rays after refraction at a distance between 5 and 10 meters. Next, we project 3D points into the second view to get the feature points. Zero mean Gaussian noise is added to feature points in both views. We vary the standard deviation of the noise from 0 to 2 pixels to evaluate the robustness of the methods. For each noise level, 50 simulations are executed.

The five-point algorithm (5pt) by Nistér (2004) that ignores refraction and the method (Ichimaru) by Ichimaru et al. (2019) are chosen as the comparison methods. Note that we carefully re-implement the method by Ichimaru et al. (2019) to make it work in general. In the following, we again name our method Virtual Epipolar Constraint – VEC. As the scale cannot be reliably recovered, we estimate the scale by $s = \sum_{i}^{N-1} \hat{d}^i/d^i$, where $d$ (and $\hat{d}$) is the Euclidean distance between two (estimated) points.

Experimental results are shown in Fig. 24. The proposed VEC algorithm does improve the camera pose initialized by the 5pt algorithm and it shows good robustness to noise. The runtime is also unsurprisingly significantly higher than that of the 5pt algorithm. Meanwhile, compared with Ichimaru et al. (2019), the proposed relative pose solver works better than theirs in terms of accuracy.

## 5.5 Scenario 2: Experiments on Real Data

We will now present results on real data sets that represents Scenario 2.

### 5.5.1 Real Data: Rubik's Cube Dataset

We first validate our algorithm on the real dataset provided by Haner and Aström (2015). The dataset was acquired with a mobile phone camera by photographing a Rubik's cube submerged in a rectangular acrylic plastic tank. By defining the global coordinate frame on the exterior facet of the tank, the authors measured manually the positions of the cube corners and the refractive plane, and marked the image correspondences by hand. To cope with the effect of the tank wall as suggested by Haner and Aström (2015), we consider a double-layer refraction model (i.e., air-plastic-water) with the refractive index of water and the plastic being 1.333 and 1.49, respectively. The proposed algorithm

<figure><img src="image_27.png" /><figcaption>**Fig. 25** Visualization of the reprojections of different solutions. Manually marked image corners are shown as white dots. The reprojections of the Haner's 5-point solver, the non-refractive solution, and the proposed solution are shown as orange crosses, magenta stars, and green plus-signs, respectively. The reprojections of the AGW, EPnP, and Ichimaru solver are omitted due to their large reprojection errors.</figcaption></figure>

<figure><img src="image_28.png" /><figcaption>**Fig. 26** Reconstructed poses: Haner's 5-point solution (in orange), no-refraction solution (in magenta), proposed solution (in green).</figcaption></figure>

can be extended to handle double-layer refraction as long as the refractive interfaces are known. In addition to the proposed method, we also include the pose estimation results solved by the 5-point algorithm (Haner and Aström, 2015) and by ignoring the refraction effects (also used as the initialization of the proposed
method), the AGW method by (Agrawal et al., 2012)<fnref n="1" />, the EPnP solution, and the Ichimaru's solution (Ichimaru et al., 2019). Their reprojections and the marked image corners overlaid on two of the images are shown in Fig. 25, and the quantitative reprojection errors are summarized in Table 1. As can be seen, the proposed method achieves the lowest reprojection errors. Similar to Haner and Åström (2015), we show the estimated camera poses in Fig. 26. As can be seen, the skewed pose estimates by the no-refraction solver can be refined to more accurate pose estimates using the proposed optimization. The result also showed that the AGW method cannot converge to good results on this dataset. This is possibly due to that they also try to estimate the normal, the distance, and the thickness of the refractive interface, which is also mentioned as a problem by Haner and Åström (2015).

### 5.5.2 Real Data: Aquarium Dataset

The aquarium Dataset is provided by Pedersen et al. (2018), where the original objective is to evaluate the precision of different camera calibration solutions. Here we use it to evaluate the accuracy of pose estimation. More specifically, we use the triangulated points based on their ray-tracing approach as 3D points, and the corresponding image coordinates as 2D features. The parameters of the refractive interface were measured by Pedersen et al. (2018). We apply our absolute pose estimation solution, the EPnP, the AGW, the 5pt and the Ichimaru algorithm to get the camera pose, and compare the obtained results with the given camera pose from the dataset. The results are shown in Table 2. Clearly, the proposed approach can estimate camera pose more accurately than the other methods, except for orientation estimation of the side-view camera where the Ichimaru solution showed slightly better performance. However, the Ichimaru solution fails for the top-view camera because their method assumes that the normal of the refractive interface is almost [0, 0, 1] as mentioned by Ichimaru et al. (2019). Therefore, for the side-view camera of which the normal of the refractive interface is [0, 0, 1], their method performs very well, while for the top-view camera that is not orthogonal to the interface, it fails completely. The 5pt algorithm gave the worst estimation for the top-view camera because it is a minimal solver that only takes 5 points in to account.²

## 5.6 Scenario 2: SfM Application

Finally, we consider three SfM applications on real data representing Scenario 2.

### 5.6.1 Real Data: Bathtub Dataset

In order to evaluate the absolute and relative pose estimation on real data, we carried out an experiment on a scene (see Fig. 27) consisting of a water filled bathtub with a submerged checkerboard and a floating checkerboard (for estimating the refractive interface and for ground truth camera poses). We used a GoPro Hero9 Black camera with narrow lens setting (f/2.5) and 5184 × 3888 pixels resolution and collected images at 4 different views.

We use the 5pt algorithm (Nistér, 2004) for initialization for the relative pose estimation, and the EPnP algorithm (Lepetit et al., 2009) for initialization for the absolute pose estimation. For comparison we estimate the reference camera pose based on the floating checkerboard using OpenCV (2022). Fig. 28(a) shows poses as camera frustum. We see that the estimated poses have moved away from the initialization and are very close to the ground truth camera poses.

Furthermore, we employed the estimated camera poses for a semi-dense 3D reconstruction of the submerged checkerboard, where the optical flow algorithm proposed by Kroeger et al. (2016) was used to establish feature correspondences. The refractive index of water is set to 1.333. The submerged checkerboard was segmented out of raw images, as demonstrated in Fig. 27. In Fig. 28(b), we show a refractive SfM reconstruction of the submerged checkerboard. The reconstruction appear smooth and as flat as expected (the bathtub bottom is not flat and the checkerboard is printed on a flexible piece of cardboard and weighted down in the water). We tried to apply a classical SfM reconstruction using the COLMAP software (Schonberger and Frahm, 2016) and AliceVision. Unfortunately, both COLMAP and AliceVision failed to give a sparse reconstruction result, whereas our method was able to reproduce the 3D structure of the submerged checkerboard. We also computed the RMSE of plane fitting using the reconstructed point cloud. As there is no comparison point cloud, to select an appropriate scale, we measured the length of square in the reconstructed point cloud and then resize the length to the real size (10 mm). The final RMSE is 0.00074.

### 5.6.2 Real Data: Stag Beetle Dataset

For evaluating the proposed absolute and relative pose estimation applied to SfM, we carried out an experiment.

<footnote n="1">¹ We extend their open-source implementation to handle the case 3 mentioned in their work
² We shuffled input points randomly. Therefore, better results might be obtained using other points for pose estimation.</footnote>
<table><caption>**Table 2** Quantitative result on Aquarium dataset (Pedersen et al., 2018): Camera pose and pose errors. Camera orientation is represented by Euler angles. Pose errors are computed by (5.1).</caption><thead><tr><th rowspan="2"></th><th rowspan="2"></th><th colspan="2">Orientation: (rad)</th><th colspan="2">Location: (cm)</th></tr><tr><th>Value</th><th>Error</th><th>Value</th><th>Error</th></tr></thead><tbody><tr><td rowspan="6">Side-view camera</td><td>PO</td><td>{-3.0638, 0.0750, 0.0912}</td><td>0.0065</td><td>{18.5471, 7.5373, -19.8116}</td><td>0.3173</td></tr><tr><td>EPnP</td><td>{3.0495, 0.0857, 0.1704}</td><td>0.1936</td><td>{18.3628, 6.4907, -16.6050}</td><td>3.0833</td></tr><tr><td>AGW</td><td>{0.0779, 0.0786, -3.0651}</td><td>0.0199</td><td>{18.6143, 7.8718, -18.0504}</td><td>1.6025</td></tr><tr><td>5pt</td><td>{0.1121, 0.0258, -3.0682}</td><td>0.0510</td><td>{17.4093, 6.8974, -20.1147}</td><td>1.3523</td></tr><tr><td>Ichimaru</td><td>{0.0922, 0.0723, -3.0625}</td><td>0.0053</td><td>{18.4796, 7.5011, -19.8930}</td><td>0.3762</td></tr><tr><td>Ground Truth</td><td>{-3.0607, 0.0737, 0.0970}</td><td></td><td>{18.5606, 7.3420, -19.5619}</td><td></td></tr><tr><td rowspan="6">Top-view camera</td><td>PO</td><td>{0.0333, 0.0160, 1.4738}</td><td>0.0067</td><td>{17.0245, 36.7280, 7.6306}</td><td>0.3248</td></tr><tr><td>EPnP</td><td>{0.0586, 0.0149, 1.6037}</td><td>0.1300</td><td>{16.0961, 33.0378, 10.7285}</td><td>5.1463</td></tr><tr><td>AGW</td><td>{0.0176, 0.0111, 1.4780}</td><td>0.0109</td><td>{17.1260, 38.5019, 7.8011}</td><td>1.4798</td></tr><tr><td>5pt</td><td>{0.3809, 0.1384, 2.9452}</td><td>1.4798</td><td>{15.5715, 3.0999, -10.6429}</td><td>38.5806</td></tr><tr><td>Ichimaru</td><td>{3.1070, -0.0177, -1.7093}</td><td>1.8073</td><td>{16.9642, 4.4794, 7.0688}</td><td>26.8317</td></tr><tr><td>Ground Truth</td><td>{0.0275, 0.0156, 1.4770}</td><td></td><td>{17.1433, 37.0295, 7.6538}</td><td></td></tr></tbody></table>

<figure><img src="image_29.png" /><figcaption>**Fig. 27** Representative raw and segmented refractive images of the bathtub dataset used for 3D reconstruction of the sub-merged checkerboard visible at the bottom of the bathtub. The floating checkerboard is used for estimation of the location of the water surface (the refractive interface) as well as for estimating ground truth camera poses.</figcaption></figure>

<figure><img src="image_30.png" /><figcaption>**Fig. 28** SfM reconstruction results of underwater checkerboard in the bathtub datasets, see images in Fig. 1 and Fig. 27. The plot (a) illustrates the initial (Init) and estimated (Est) camera poses together with reference (Ref) camera poses estimated from the floating checkerboard - the ground truth camera poses. In (b), the left (top view) and right image (side view) shows SfM reconstruction based on refractive SfM using our method for Scenario 2.</figcaption></figure>
<figure><img src="image_31.png" /><figcaption>**Fig. 29** Representative raw and segmented refractive images of the real stag beetle dataset used for 3D reconstruction.</figcaption></figure>

iment where we used a NIKON D40X (focal length: 34mm, resolution: 3872 × 2592 pixels) to photograph a stag beetle enclosed in crystal resin, as shown in Fig. 29, and collected images in 30 different views. The refractive index of the crystal resin is set to 1.6. We segmented the refractive parts from raw images, as demonstrated in Fig. 29, on which the proposed algorithms were then executed. Four ArUco (Garrido-Jurado et al., 2016) fiducial markers were glued on the front surface of the resin to estimate the refractive interface in the world coordinate frame (the first view defines the world frame). We initialize the camera poses from the background of the scene using AliceVision (AliceVision, 2022; Moulon et al., 2012).

For comparison, we compute a classical SfM reconstruction using the COLMAP software (Schönberger and Frahm, 2016). In Fig. 30, we show two classical SfM reconstructions and a refractive SfM reconstruction of the stag beetle from Fig. 29. Notice that the classical reconstruction fails to get the depth correct and produces an almost flat reconstruction, whereas our method more faithfully reproduce the 3D structure of the stag beetle (especially the legs and other attachments on the body), see Fig. 1 for a side-view of the real stag beetle. However, the refractive reconstruction is slightly more noisy.

### 5.6.3 Whiskey Bottle Dataset

We carried out an additional experiment where we used a mobile phone camera to photograph a cylindrical whiskey bottle container placed at the bottom of a bathtub

filled with water, as shown in Fig. 31. Similarly, we segmented the refractive parts from raw images on which the proposed algorithms were then executed. In Fig. 32, we show the submerged whiskey bottle and tray reconstructed using our incremental RSfM, AliceVision, and COLMAP. It is observed from Fig. 32 that AliceVision and COLMAP performs well on this dataset. Visually speaking, AliceVision produce even better results than ours. We consider the possible reasons could be that the camera intrinsic parameters are also optimized during SfM. It has been shown in previous work (Lavest et al., 2000; Łuczyński et al., 2017) that the effects caused by refraction can be reduced to some extent by modeling them as additional radial distortions and adjusting the radial distortion parameters. This compensation may work well especially when the camera does not move a lot in the scene. Nonetheless, our RSfM also give promising reconstruction results, which validates our pose estimation solution.

## 6 Conclusion

In this paper, we have defined the theoretical framework for absolute and relative camera pose estimation under refraction given the refractive interface parameters and 3D/2D-2D feature correspondences. According to the coordinate frame where the refractive interface is given, we define two scenarios where the first scenario corresponds to the underwater scenarios and the second one represents the scenario when an object is photographed through a refractive surface such as amber
<figure><img src="image_32.png" /><figcaption>**Fig. 30** SfM reconstructions on the stag beetle dataset: (Left) SfM reconstruction based on refractive SfM using our method for Scenario 2. (Middle) SfM reconstruction based on perspective camera using *AliceVision* (2022). (Right) SfM reconstruction based on perspective camera using *COLMAP* (Schonberger and Frahm, 2016).</figcaption></figure>

<figure><img src="image_33.png" /><figcaption>**Fig. 31** Representative raw and segmented refractive images of the real whisky bottle dataset used for 3D reconstruction.</figcaption></figure>

or an aquarium. For each scenario, we have derived its absolute and relative pose estimation solutions which are primarily built upon our proposed iterative refinements. More specially, for the underwater scenario, we can also derive non-iterative solutions that are analogue to classical pose estimation methods (i.e. PnP or GEC) in non-refractive space, which can work in low-noise situations. We have demonstrated the superiority of our solution compared to State-of-the-art on synthetic and

real data. We also show that, by integrating our pose estimation solutions in a RSfM framework, the RSfM can reconstruct objects under refraction properly and show better results compared with classical SfM methods on selected dataset.

For scenario 1, the underwater scenario, we show in *Section 5.2.1* that our approach provides superior results in camera pose estimation in real applications. Similarly, in *Section 5.3*, we study the effect of refrac-
<figure><img src="image_34.png" /><figcaption>Fig. 32 Representative side-view images of the real whisky bottle dataset: (Left) the reconstruction results using our incremental RSfM; (Middle) the reconstruction results using AliceVision (2022); (Right) the reconstruction results using COLMAP (Schonberger and Frahm, 2016).</figcaption></figure>

tion on SfM reconstruction in the underwater scenario and confirm previous findings (Jordt and Koch, 2013; Chadebecq et al., 2019).

For scenario 2, we show in Section 5.5 the superiority of our approach to camera pose estimation in real applications. In Section 5.6 we show that using RSfM on real applications provides improved 3D reconstructions of the scenes.

As for the future work, we will consider the following points:

- For scenario 2, we are still missing proper initialization algorithms. We will consider to propose initialization algorithms targeting this issue.

- We will consider to integrate our solutions in a robust estimation framework so that outliers can be handled.

- We will integrate our solutions into some mature SfM framework like AliceVision so that its efficient and robust pipeline can be used to ameliorate the quality and speed of RSfM.

- As our project focuses on reconstructing insects inside Amber, we will also consider to extend our solutions to handle non-planar refractive interfaces, e.g. (Xiong and Heidrich, 2021b).

- Finally, as we point out in the introduction, feature matching under refraction is also a challenging task, since classical geometric verification models are invalid. To the best of our knowledge, there is not yet a well-established solution to identify erroneous matches under refraction. Future work will consider how to handle the erroneous matches in feature matching.

## A Jacobian Matrices related to Pose Estimation under Refraction

This appendix section gives the Jacobian matrices derived for pose estimation under refraction.

### A.1 Scenario 1

#### A.1.1 Relative Pose Refinement

The residual of the epipolar cost and its analytical Jacobian are given as follows:

$$ 
\begin{align} 
\mathcal{E} &= \mathbf{u}'^{\top} \mathbf{R}^{\top} (\mathbf{e})^{\top} \mathbf{u}, \nonumber \\
\mathcal{J} \mathcal{E}_{\mathbf{x}} &= [\mathcal{J} \mathcal{E}_{\mathbf{R}}, \mathcal{J} \mathcal{E}_{\mathbf{c}}, \mathbf{0}], \nonumber \\
\mathcal{J} \mathcal{E}_{\mathbf{R}} &= \mathbf{u}'^{\top} \mathbf{R}^{\top} \mathbf{u}^{\top} \mathbf{R} \mathcal{c}_{\mathcal{H}_2} - \mathbf{u}'^{\top} (\mathbf{R}^{\top} \mathbf{u}^{\top} \mathbf{e})^{\top}, \tag{A.1} \\
\mathcal{J} \mathcal{E}_{\mathbf{c}} &= -\mathbf{u}'^{\top} \mathbf{R}^{\top} \mathbf{u}^{\top}, \nonumber
\end{align} 
$$

where $\mathbf{e} = \mathbf{R} \mathbf{c}_{V_2} + \mathbf{c} - \mathbf{c}_{V_1}$, $\mathbf{x} = [\boldsymbol{\xi}^{\top}, \mathbf{v}^{\top}, \mathbf{w}^{i^{\top}}]^{\top}$, $i = 1, \dots, N$.

Next, we derive the residual function of the reprojection cost and the corresponding Jacobian matrix as follows:

$$ 
\begin{align} 
\mathcal{R} &= \begin{bmatrix} f_V \frac{\mathbf{p}_{V,1}}{\mathbf{p}_{V,3}} - f_V \frac{\mathbf{r}_{V,1}}{\mathbf{r}_{V,3}} \\ f_V \frac{\mathbf{p}_{V,2}}{\mathbf{p}_{V,3}} - f_V \frac{\mathbf{r}_{V,2}}{\mathbf{r}_{V,3}} \end{bmatrix}, \nonumber \\
\mathcal{J} \mathcal{R}_{\mathbf{x}} &= \mathbf{J}_{\mathbf{p}_V} \mathbf{J}_{\mathbf{x}}^{\mathrm{v}}, \nonumber \\
\mathbf{J}_{\mathbf{p}_V} &= \begin{bmatrix} \frac{f_V}{\mathbf{p}_{V,3}} & 0 & -f_V \frac{\mathbf{p}_{V,1}}{\mathbf{p}_{V,3}} \\ 0 & \frac{f_V}{\mathbf{p}_{V,3}} & -f_V \frac{\mathbf{p}_{V,2}}{\mathbf{p}_{V,3}} \end{bmatrix}, \tag{A.2} \\
\mathbf{J}_{\mathbf{x}}^{\mathrm{v}} &= [\mathbf{R}_{\mathbf{V}}^{\top} (\mathbf{R}^{\top} (\mathbf{p} - \mathbf{c}))^{\top} - \mathbf{R}_{\mathbf{V}}^{\top} \mathbf{R}^{\top} \mathbf{R}_{\mathbf{V}}^{\top} \mathbf{R}^{\top}] , \nonumber
\end{align} 
$$

where $\mathbf{p}_V = \mathbf{R}_{\mathbf{V}}^{\top} (\mathbf{R}^{\top} (\mathbf{p} - \mathbf{c}) - \mathbf{c}_{V_2})$, $\mathbf{r}_V = \mathbf{R}_{\mathbf{V}}^{\top} \mathbf{r}$, $\mathbf{r}$ denotes the refractive ray of feature $\mathbf{u}$, $f_V$ denotes the focal length of the virtual camera, $\mathbf{R}_V$ and $\mathbf{c}_V$ are the rotation and location of the virtual camera coordinate frame with respect to the camera coordinate frame.
## A.2 Scenario 2

### A.2.1 Absolute Pose Refinement

By using the chain rule, we have the Jacobian matrix w.r.t **R** and **C** as

$$ 
\begin{aligned}
\mathbf{J}\mathcal{L}_{\mathbf{R}} &= \mathbf{J}\mathbf{I}_{\mathbf{W}_1}\mathbf{J}_{\mathbf{R}}^{\mathbf{I}_{\mathbf{W}_1}} + \mathbf{J}_{\mathbf{W}_2}\mathbf{J}_{\mathbf{R}}^{\mathbf{I}_{\mathbf{W}_2}}, \\
\mathbf{J}_{\mathbf{c}} &= \mathbf{J}_{\mathbf{I}_{\mathbf{W}_1}}\mathbf{J}_{\mathbf{c}}^{\mathbf{I}_{\mathbf{W}_1}} + \mathbf{J}_{\mathbf{I}_{\mathbf{W}_2}}\mathbf{J}_{\mathbf{c}}^{\mathbf{I}_{\mathbf{W}_2}}.
\end{aligned}
\quad (A.3) $$

where $\mathbf{J}_{\mathbf{I}_{\mathbf{W}_1}} = -\mathbf{p}_{\mathbf{W}}^{\wedge}$ and $\mathbf{J}_{\mathbf{I}_{\mathbf{W}_2}} = \mathbf{I}$. Since **R** and **c** are embedded in $\mathbf{I}_{\mathbf{W}}$, again, by the chain rule, we have:

$$ 
\begin{aligned}
\mathbf{J}_{\mathbf{R}}^{\mathbf{I}_{\mathbf{W}_1}} &= \mathbf{J}_{\mathbf{q}_W}^{\mathbf{I}_{\mathbf{W}_1}} \mathbf{J}_{\mathbf{I}_W}^{\mathbf{q}_W} \mathbf{J}_{\mathbf{R}}^{\mathbf{I}_W} + \mathbf{J}_{\mathbf{r}_W}^{\mathbf{I}_{\mathbf{W}_1}} \mathbf{J}_{\mathbf{I}_W}^{\mathbf{r}_W} \mathbf{J}_{\mathbf{R}}^{\mathbf{I}_W}, \\
\mathbf{J}_{\mathbf{c}}^{\mathbf{I}_{\mathbf{W}_1}} &= \mathbf{J}_{\mathbf{q}_W}^{\mathbf{I}_{\mathbf{W}_1}} \mathbf{J}_{\mathbf{w}_{c_C}}^{\mathbf{q}_W} \\
\mathbf{J}_{\mathbf{R}}^{\mathbf{I}_{\mathbf{W}_2}} &= \mathbf{J}_{\mathbf{q}_W}^{\mathbf{I}_{\mathbf{W}_2}} \mathbf{J}_{\mathbf{I}_W}^{\mathbf{q}_W} \mathbf{J}_{\mathbf{R}}^{\mathbf{I}_W} + \mathbf{J}_{\mathbf{r}_W}^{\mathbf{I}_{\mathbf{W}_2}} \mathbf{J}_{\mathbf{I}_W}^{\mathbf{r}_W} \mathbf{J}_{\mathbf{R}}^{\mathbf{I}_W}, \\
\mathbf{J}_{\mathbf{c}}^{\mathbf{I}_{\mathbf{W}_2}} &= \mathbf{J}_{\mathbf{q}_W}^{\mathbf{I}_{\mathbf{W}_2}} \mathbf{J}_{\mathbf{w}_{c_C}}^{\mathbf{q}_W} \\
\mathbf{J}_{\mathbf{I}_W}^{\mathbf{q}_W} &= \begin{bmatrix} -\mathbf{W}\mathbf{c}_{C,3} & 0 & \mathbf{W}\mathbf{c}_{C,1} \\ 0 & -\mathbf{W}\mathbf{c}_{C,3} & \mathbf{W}\mathbf{c}_{C,2} \\ 0 & 0 & 0 \\ 0 & 0 & 1 \end{bmatrix},
\end{aligned}
\quad (A.4) $$

where

$$ 
\begin{aligned}
\mathbf{J}_{\mathbf{q}_W}^{\mathbf{I}_{\mathbf{W}_1}} &= [\mathbf{0} \quad \mathbf{r}_W], \quad \mathbf{J}_{\mathbf{r}_W}^{\mathbf{I}_{\mathbf{W}_1}} = [\mathbf{q}_{W,4}\mathbf{I} \quad \mathbf{0}], \\
\mathbf{J}_{\mathbf{q}_W}^{\mathbf{I}_{\mathbf{W}_2}} &= [-\mathbf{r}_{W}^{\wedge} \quad \mathbf{0}], \quad \mathbf{J}_{\mathbf{r}_W}^{\mathbf{I}_{\mathbf{W}_2}} = [\mathbf{q}_{W}^{\wedge} \quad \mathbf{0}], \\
\mathbf{J}_{\mathbf{I}_W}^{\mathbf{q}_W} &= \begin{bmatrix} -\mathbf{W}\mathbf{c}_{C,3} & 0 & \mathbf{W}\mathbf{c}_{C,1} \\ 0 & -\mathbf{W}\mathbf{c}_{C,3} & \mathbf{W}\mathbf{c}_{C,2} \\ 0 & 0 & 0 \\ 0 & 0 & 1 \end{bmatrix},
\end{aligned}
$$

$$ 
\begin{aligned}
\mathbf{J}_{\mathbf{w}_{c_C}}^{\mathbf{q}_W} &= \begin{bmatrix} \mathbf{i}_{W,3} & 0 & -\mathbf{i}_{W,1} \\ 0 & \mathbf{i}_{W,3} & -\mathbf{i}_{W,2} \\ 0 & 0 & 0 \\ 0 & 0 & 0 \end{bmatrix}, \\
\mathbf{J}_{\mathbf{I}_W}^{\mathbf{r}_W} &= \begin{bmatrix} \lambda & 0 & 0 \\ 0 & \lambda & 0 \\ 0 & 0 & \frac{\lambda^2 \mathbf{i}_{W,3}}{\mathbf{r}_{W,3}} \\ 0 & 0 & 0 \end{bmatrix}, \\
\mathbf{J}_{\mathbf{R}}^{\mathbf{I}_W} &= (\mathbf{i}_W)^{\wedge}.
\end{aligned}
\quad (A.5) $$

### A.2.2 Relative Pose Refinement

Recall the refractive epipolar constraint given as

$$ \mathcal{E} = \mathbf{r}_{W}^{2}^{\top} (\mathbf{W}^{2}\mathbf{c}_{V}^{2} - \mathbf{W}^{2}\mathbf{c}_{V}^{1})^{\wedge} \mathbf{r}_{W}^{1}, \quad (A.6) $$

using the chain rule, we have the Jacobian matrices w.r.t the variation of **R**, **c**, and **p** as

$$ 
\begin{aligned}
\mathbf{J}\mathcal{E}_{\mathbf{x}} &= [\mathbf{J}\mathcal{E}_{\mathbf{R}}, \mathbf{J}\mathcal{E}_{\mathbf{c}}, \mathbf{0}], \\
\mathbf{J}\mathcal{E}_{\mathbf{R}} &= \mathbf{J}\mathcal{E}_{\mathbf{r}_W}^{2} \mathbf{J}_{\mathbf{R}}^{\mathbf{r}_W} + \mathbf{J}\mathcal{E}_{\mathbf{w}_{c_V}^{2}} \mathbf{J}_{\mathbf{R}}^{\mathbf{w}_{c_V}^{2}}, \\
\mathbf{J}\mathcal{E}_{\mathbf{c}} &= \mathbf{J}\mathcal{E}_{\mathbf{w}_{c_V}^{2}} \mathbf{J}_{\mathbf{c}}^{\mathbf{w}_{c_V}^{2}}.
\end{aligned}
\quad (A.7) $$

where

$$ 
\begin{aligned}
\mathbf{J}\mathcal{E}_{\mathbf{r}_W}^{2} &= ((\mathbf{W}^{2}\mathbf{c}_{V}^{2} - \mathbf{W}^{2}\mathbf{c}_{V}^{1})^{\top} \mathbf{r}_{W}^{1})^{\top}, \\
\mathbf{J}\mathcal{E}_{\mathbf{w}_{c_V}^{2}} &= -\mathbf{r}_{W}^{2}^{\top} (\mathbf{r}_{W}^{1})^{\wedge}, \\
\mathbf{J}_{\mathbf{R}}^{\mathbf{r}_W} &= \mathbf{J}_{\mathbf{I}_W}^{\mathbf{r}_W} \mathbf{J}_{\mathbf{R}}^{\mathbf{I}_W}, \\
\mathbf{J}_{\mathbf{R}}^{\mathbf{w}_{c_V}^{2}} &= \mathbf{J}_{\mathbf{I}_W}^{\mathbf{w}_{c_V}^{2}} \mathbf{J}_{\mathbf{R}}^{\mathbf{I}_W}, \\
\mathbf{J}_{\mathbf{c}}^{\mathbf{w}_{c_V}^{2}} &= \mathbf{J}_{\mathbf{w}_{c_C}}^{\mathbf{w}_{c_V}^{2}}.
\end{aligned}
\quad (A.8) $$

Similarly, for the reprojection error given as

$$ \mathcal{R} = \begin{bmatrix} fV \frac{\mathbf{p}_{\mathbf{V},1}}{\mathbf{p}_{\mathbf{V},3}} - fV \frac{\mathbf{r}_{\mathbf{V},1}}{\mathbf{r}_{\mathbf{V},3}} \\ fV \frac{\mathbf{p}_{\mathbf{V},2}}{\mathbf{p}_{\mathbf{V},3}} - fV \frac{\mathbf{r}_{\mathbf{V},2}}{\mathbf{r}_{\mathbf{V},3}} \end{bmatrix}, \quad (A.9) $$

using the chain rule, we have the Jacobian matrices w.r.t **R**, **c**, and **p** as

$$ 
\begin{aligned}
\mathbf{J}\mathcal{R}_{\mathbf{x}} &= [\mathbf{J}\mathcal{R}_{\mathbf{R}}, \mathbf{J}\mathcal{R}_{\mathbf{c}}, \mathbf{J}\mathcal{R}_{\mathbf{p}}], \\
\mathbf{J}\mathcal{R}_{\mathbf{R}} &= \mathbf{J}\mathcal{R}_{\mathbf{p}_V} \mathbf{J}_{\mathbf{R}}^{\mathbf{p}_V} + \mathbf{J}\mathcal{R}_{\mathbf{r}_V} \mathbf{J}_{\mathbf{R}}^{\mathbf{r}_V}, \\
\mathbf{J}\mathcal{R}_{\mathbf{c}} &= \mathbf{J}\mathcal{R}_{\mathbf{p}_V} \mathbf{J}_{\mathbf{c}}^{\mathbf{p}_V} + \mathbf{J}\mathcal{R}_{\mathbf{r}_V} \mathbf{J}_{\mathbf{c}}^{\mathbf{r}_V}, \\
\mathbf{J}\mathcal{R}_{\mathbf{p}} &= \mathbf{J}\mathcal{R}_{\mathbf{p}_V} \mathbf{J}_{\mathbf{p}}^{\mathbf{p}_V},
\end{aligned}
\quad (A.10) $$

where

$$ 
\begin{aligned}
\mathbf{J}\mathcal{R}_{\mathbf{p}_V} &= \begin{bmatrix} fV \frac{1}{\mathbf{p}_{\mathbf{V},3}} & 0 & -fV \frac{\mathbf{p}_{\mathbf{V},1}}{\mathbf{p}_{\mathbf{V},3}} \\ 0 & fV \frac{1}{\mathbf{p}_{\mathbf{V},3}} & -fV \frac{\mathbf{r}_{\mathbf{V},1}}{\mathbf{r}_{\mathbf{V},3}} \end{bmatrix}, \\
\mathbf{J}\mathcal{R}_{\mathbf{r}_V} &= \begin{bmatrix} -fV \frac{1}{\mathbf{r}_{\mathbf{V},3}} & 0 & fV \frac{\mathbf{r}_{\mathbf{V},1}}{\mathbf{r}_{\mathbf{V},3}} \\ 0 & -fV \frac{1}{\mathbf{r}_{\mathbf{V},3}} & fV \frac{\mathbf{r}_{\mathbf{V},1}}{\mathbf{r}_{\mathbf{V},3}} \end{bmatrix}, \\
\mathbf{J}_{\mathbf{R}}^{\mathbf{p}_V} &= \mathbf{J}_{\mathbf{w}_{c_V}^{2}} \mathbf{J}_{\mathbf{R}}^{\mathbf{w}_{c_V}^{2}}, \quad \mathbf{J}_{\mathbf{R}}^{\mathbf{r}_V} = \mathbf{J}_{\mathbf{q}_W}^{\mathbf{r}_W} \mathbf{J}_{\mathbf{I}_W}^{\mathbf{q}_W} \mathbf{J}_{\mathbf{R}}^{\mathbf{I}_W}, \quad \mathbf{J}_{\mathbf{p}}^{\mathbf{p}_V} = \mathbf{I}, \\
\mathbf{J}_{\mathbf{c}}^{\mathbf{p}_V} &= \mathbf{J}_{\mathbf{w}_{c_V}^{2}} \mathbf{J}_{\mathbf{w}_{c_C}}^{\mathbf{w}_{c_V}^{2}}, \quad \mathbf{J}_{\mathbf{c}}^{\mathbf{r}_V} = \mathbf{J}_{\mathbf{q}_W}^{\mathbf{r}_W} \mathbf{J}_{\mathbf{w}_{c_C}}^{\mathbf{w}_{c_V}^{2}}.
\end{aligned}
\quad (A.11) $$

Some intermediate Jacobian matrices in (A.8) and (A.11) are given as follows:

$$ 
\begin{aligned}
\mathbf{J}_{\mathbf{R}}^{\mathbf{r}_W} &= \mathbf{J}_{\mathbf{I}_W}^{\mathbf{r}_W} \mathbf{J}_{\mathbf{R}}^{\mathbf{I}_W}, \\
\mathbf{J}_{\mathbf{R}}^{\mathbf{w}_{c_V}^{2}} &= \mathbf{J}_{\mathbf{I}_W}^{\mathbf{w}_{c_V}^{2}} \mathbf{J}_{\mathbf{R}}^{\mathbf{w}_{c_V}^{2}}, \\
\mathbf{J}_{\mathbf{c}}^{\mathbf{w}_{c_V}^{2}} &= \mathbf{J}_{\mathbf{w}_{c_C}}^{\mathbf{w}_{c_V}^{2}},
\end{aligned}
\quad (A.12) $$

where

$$ 
\begin{aligned}
\mathbf{J}_{\mathbf{w}_{c_C}}^{\mathbf{q}_W} &= \begin{bmatrix} 1 & 0 & -\frac{\mathbf{i}_{W,1}}{\mathbf{i}_{W,3}} \\ 0 & 1 & -\frac{\mathbf{i}_{W,2}}{\mathbf{i}_{W,3}} \\ 0 & 0 & 0 \end{bmatrix}, \\
\mathbf{J}_{\mathbf{I}_W}^{\mathbf{q}_W} &= \begin{bmatrix} -\frac{\mathbf{w}\mathbf{c}_{C,3}}{\mathbf{i}_{W,3}} & 0 & \mathbf{W}\mathbf{c}_{C,3} & \frac{\mathbf{i}_{W,1}}{\mathbf{i}_{W,3}} \\ 0 & -\frac{\mathbf{w}\mathbf{c}_{C,3}}{\mathbf{i}_{W,3}} & \mathbf{W}\mathbf{c}_{C,3} & \frac{\mathbf{i}_{W,1}}{\mathbf{i}_{W,3}} \\ 0 & 0 & 0 & 0 \end{bmatrix}, \\
\mathbf{J}_{\mathbf{I}_W}^{\mathbf{r}_W} &= \begin{bmatrix} \lambda & 0 & 0 \\ 0 & \lambda & 0 \\ 0 & 0 & \frac{\lambda^2 \mathbf{i}_{W,3}}{\gamma} \end{bmatrix}, \\
\mathbf{J}_{\mathbf{w}_{c_V}^{2}}^{\mathbf{c}^i V} &= \begin{bmatrix} 0 & 0 & 0 \\ 0 & 0 & 0 \\ 0 & 0 & \frac{-\mathbf{w}\mathbf{c}_{C,3}(1-\lambda^2)}{\lambda \mathbf{i}_{W,3}^2 \gamma} \end{bmatrix}, \\
\mathbf{J}_{\mathbf{w}_{c_C}^{2}}^{\mathbf{c}^i V} &= \begin{bmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & \frac{1}{\lambda \mathbf{i}_{W,3}} \gamma \end{bmatrix}, \\
\mathbf{J}_{\mathbf{R}}^{\mathbf{I}_W} &= (\mathbf{q}_W)^{\wedge}, \quad \mathbf{J}_{\mathbf{w}_{c_V}^{2}}^{\mathbf{p}_V} = -\mathbf{I}, \quad \mathbf{J}_{\mathbf{q}_W}^{\mathbf{r}_W} = \mathbf{I}.
\end{aligned}
\quad (A.13) $$

where $\gamma = \sqrt{(1 - \lambda^2) + \lambda^2 \mathbf{i}_{W,3}^2}$.
# B Additional Experiments

## B.1 Setup of Synthetic Simulation

Fig. 33 show how the data are generated to evaluate the absolute pose solvers in the non planar and planar cases. Fig. 34 show how the data are generated to evaluate the relative pose solvers.

<figure><img src="image_35.png" /></figure>

**Fig. 33** Simulation setup for absolute pose estimation: The left figure shows the setup for generating data (world points, image points, camera pose) under refraction for the non planar case. The right figure shows the setup for generating data for the planar case.

<figure><img src="image_36.png" /></figure>

**Fig. 34** Simulation setup for relative pose estimation: The figure shows the setup for generating data (world points, image points, camera pose) under refraction.

## B.2 Scenario 1: Absolute Pose Solver

As mentioned in Section 5.1.1, we found that

- UPnP Kneip et al. (2014) unfortunately does not perform well even in non planar case without noise perturbation, which is shown in Fig. 35.

- gPnP Kneip et al. (2013), gP+s Ventura et al. (2014), and MMPPE Miraldo et al. (2018) unfortunately does not perform well even in planar case without noise perturbation, which is shown in Fig. 36.

<figure><img src="image_37.png" /></figure>

**Fig. 35** Comparison of absolute pose solvers in non planar case. Note that we set the standard deviation of the Gaussian noise to 0.

<figure><img src="image_38.png" /></figure>

**Fig. 36** Comparison of absolute pose solvers in planar case. Note that we set the standard deviation of the Gaussian noise to 0.

## B.3 Scenario 1: Relative Pose Solver

As mentioned in Section 5.1.2, we intended to add other GEC (Generalized Epipolar Constraint Pless (2003)) such as ge Kneip and Li (2014) and 6pt Stewenius et al. (2005). However, ge and 6pt do not show satisfactory performance in noise-free condition, as shown in Fig. 37.

<figure><img src="image_39.png" /></figure>

**Fig. 37** Comparison of GEC solvers in noise-free case. Note that 17pt gives very accurate results so that its box is very tiny in the box plots.

# References

- Absil PA, Mahony R, Sepulchre R (2009) Optimization algorithms on matrix manifolds. Princeton University Press
Agrawal A, Ramalingam S, Taguchi Y, Chari V (2012) A theory of multi-layer flat refractive geometry. In: Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, IEEE, pp 3346-3353

AliceVision (2022) Photogrammetric computer vision framework. https://alicevision.org/

Arun KS, Huang TS, Blostein SD (1987) Least-squares fitting of two 3-d point sets. IEEE Transactions on Pattern Analysis and Machine Intelligence PAMI-9(5):698-700, DOI 10.1109/TPAMI.1987.4767965

Cassidy M, Mélou J, Quéau Y, Lauze F, Durou JD (2020) Refractive multi-view stereo. In: International Conference on 3D Vision (3DV 2020)

Chadebecq F, Vasconcelos F, Dwyer G, Lacher R, Ourselin S, Vercauteren T, Stoyanov D (2017) Refractive structure-from-motion through a flat refractive interface. In: Proceedings of the IEEE International Conference on Computer Vision, pp 5315-5323

Chadebecq F, Vasconcelos F, Lacher R, Maneas E, Desjardins A, Ourselin S, Vercauteren T, Stoyanov D (2019) Refractive two-view reconstruction for underwater 3d vision. International Journal of Computer Vision DOI 10.1007/s11263-019-01218-9

Chang YJ, Chen T (2011) Multi-view 3d reconstruction for scenes under the refractive plane with known vertical direction. In: 2011 International Conference on Computer Vision, IEEE, pp 351-358

Chari V, Sturm P (2009) Multiple-view geometry of the refractive plane. In: BMVC 2009-20th British Machine Vision Conference, The British Machine Vision Association (BMVA), pp 1-11

Ferraz L, Binefa X, Moreno-Noguer F (2014) Very fast solution to the pnp problem with algebraic outlier rejection. In: Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, pp 501-508

Fischler MA, Bolles RC (1981) Random sample consensus: a paradigm for model fitting with applications to image analysis and automated cartography. Communications of the ACM 24(6):381-395

Fragoso V, DeGol J, Hua G (2020) gdls*: Generalized pose-and-scale estimation given scale and gravity priors. In: Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pp 2210-2219

Gao XS, Hou XR, Tang J, Cheng HF (2003) Complete solution classification for the perspective-three-point problem. IEEE transactions on pattern analysis and machine intelligence 25(8):930-943

Garrido-Jurado S, Munoz-Salinas R, Madrid-Cuevas FJ, Medina-Carnicer R (2016) Generation of fiducial marker dictionaries using mixed integer linear programming. Pattern Recognition 51:481-491

gP+s (2014) https://github.com/jonathanventura/genposeandscale

Grossberg MD, Nayar SK (2001) A general imaging model and a method for finding its parameters. In: Proceedings of the IEEE International Conference on Computer Vision. ICCV 2001, IEEE, vol 2, pp 108-115

Hadfield S, Lebeda K, Bowden R (2018) Hard-pnp: Pnp optimization using a hybrid approximate representation. IEEE transactions on pattern analysis and machine intelligence 41(3):768-774

Haner S, Åström K (2015) Absolute pose for cameras under flat refractive interfaces. In: Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, pp 1428-1436

Hartley RI, Zisserman A (2004) Multiple View Geometry in Computer Vision, 2nd edn. Cambridge University Press, ISBN: 0521540518

Hedrick BP, Heberling JM, Meineke EK, Turner KG, Gassa CJ, Park DS, Kennedy J, Clarke JA, Cook JA, Blackburn DC, Edwards SV, Davis CC (2020) Digitization and the future of natural history collections. BioScience 70(3):243-251

Hesch JA, Roumeliotis SI (2011) A direct least-squares (dls) method for pnp. In: Proceedings of the IEEE International Conference on Computer Vision, IEEE, pp 383-390

Hu X, Lauze F, Pedersen KS, Mélou J (2021) Absolute and relative pose estimation in refractive multi view. In: Proceedings of the IEEE/CVF International Conference on Computer Vision, pp 2569-2578

Hu X, Lauze F, Pedersen KS (2022a) RefractiveSfM. https://github.com/diku-dk/RefractiveSfM

Hu X, Lauze F, Pedersen KS, Quéau Y (2022b) DIKU Refractive Scenes Dataset 2022. Data DOI: http://doi.org/10.17894/ucph.5d1b9bea-b105-4d43-aefb-c53df7806c2a

Ichimaru K, Taguchi Y, Kawasaki H (2019) Unified underwater structure-from-motion. In: Proceedings - 2019 International Conference on 3D Vision, 3DV 2019, pp 524-532

Jordt A (2013) Underwater 3d reconstruction based on physical models for refraction and underwater light propagation. PhD thesis, URL https://macau.uni-kiel.de/receive/diss_mods_00014162

Jordt A, Koch R (2012) Refractive calibration of underwater cameras. In: Proceedings of the European Conference on Computer Vision, Springer, pp 846-859

Jordt A, Koch R (2013) Refractive structure-from-motion on underwater images. In: Proceedings of the IEEE International Conference on Computer Vision, pp 57-64

Jordt A, Köser K, Koch R (2016) Refractive 3d reconstruction on underwater images. Methods in Oceanography 15:90-113

Kang L, Wu L, Yang YH (2012a) Experimental study of the influence of refraction on underwater three-dimensional reconstruction using the svp camera model. Applied optics 51(31):7591-7603

Kang L, Wu L, Yang YH (2012b) Two-view underwater structure and motion for cameras under flat refractive interfaces. In: Proceedings of the European Conference on Computer Vision, Springer, pp 303-316

Kneip L, Furgale P (2014) Opengv: A unified and generalized approach to real-time calibrated geometric vision. In: 2014 IEEE International Conference on Robotics and Automation (ICRA), IEEE, pp 1-8

Kneip L, Li H (2014) Efficient computation of relative pose for multi-camera systems. In: Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, pp 446-453

Kneip L, Scaramuzza D, Siegwart R (2011) A novel parametrization of the perspective-three-point problem for a direct computation of absolute camera position and orientation. In: Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, pp 2969-2976

Kneip L, Furgale P, Siegwart R (2013) Using multi-camera systems in robotics: Efficient solutions to the npnp problem. In: Proceedings of the IEEE International Conference on Robotics and Automation, IEEE, pp 3770-3776

Kneip L, Li H, Seo Y (2014) Upnp: An optimal o (n) solution to the absolute pose problem with universal applicability. In: Proceedings of the European Conference on Computer Vision, Springer, pp 127-142
Kroeger T, Timofte R, Dai D, Van Gool L (2016) Fast optical flow using dense inverse search. In: Proceedings of the European Conference on Computer Vision, Springer, pp 471–488

Kukelova Z, Bujnak M, Pajdla T (2008) Automatic generator of minimal problem solvers. In: Proceedings of the European Conference on Computer Vision, Springer, pp 302–315

Lavest JM, Rives G, Lapresté JT (2000) Underwater camera calibration. In: Proceedings of the European Conference on Computer Vision, Springer, pp 654–668

Lee GH, Li B, Pollefeys M, Fraundorfer F (2016) Minimal solutions for pose estimation of a multi-camera system. In: 16th International Symposium of Robotics Research, ISRR 2013, Springer Verlag, pp 521–538

Lepetit V, Moreno-Noguer F, Fua P (2009) Epnp: An accurate o (n) solution to the pnp problem. International journal of computer vision 81(2):155

Li H, Hartley R, Kim Jh (2008) A linear approach to motion estimation using generalized camera models. In: Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, IEEE, pp 1–8

Lowe DG (2004) Distinctive image features from scale-invariant keypoints. International journal of computer vision 60(2):91–110

Łuczyński T, Pfingsthorn M, Birk A (2017) The pinax-model for accurate and efficient refraction correction of underwater cameras in flat-pane housings. Ocean Engineering 133:9–22

Ma Y, Soatto S, Košecká J, Sastry S (2004) An invitation to 3-D vision. From Images to Geometric Models. Interdisciplinary Applied Mathematics, Springer

Miraldo P, Dias T, Ramalingam S (2018) A minimal closed-form solution for multi-perspective pose estimation using points and lines. In: Proceedings of the European Conference on Computer Vision, pp 474–490

MMPPE (2018) https://github.com/pmiraldo/MinimalMultiPerspectivePose

Moulon P, Monasse P, Marlet R (2012) Adaptive structure from motion with a contrario model estimation. In: Proceedings of the Asian Computer Vision Conference, Springer Berlin Heidelberg, pp 257–270, DOI 10.1007/978-3-642-37447-0-20

Mouragnon E, Lhuillier M, Dhome M, Dekeyser F, Sayd P (2009) Generic and real-time structure from motion using local bundle adjustment. Image and Vision Computing 27(8):1178–1193

NHMD (2022a) Amber collection, Natural History Museum of Denmark (NHMD). https://samlinger.snm.ku.dk/en/dry-and-wet-collections/zoology/entomology/amber-collection/

NHMD (2022b) Herpetology collection, Natural History Museum of Denmark (NHMD). https://samlinger.snm.ku.dk/en/dry-and-wet-collections/zoology/herpetology-collection/

Nistér D (2004) An efficient solution to the five-point relative pose problem. IEEE transactions on pattern analysis and machine intelligence 26(6):756–770

Nistér D, Stewénius H (2007) A minimal solution to the generalised 3-point pose problem. Journal of Mathematical Imaging and Vision 27(1):67–79

Oleari F, Kallasi F, Rizzini DL, Aleotti J, Caselli S (2015) An underwater stereo vision system: From design to deployment and dataset acquisition. In: OCEANS 2015-Genova, IEEE, pp 1–6

OpenCV (2022) Open source computer vision library. https:// opengl.org/

Pedersen M, Bengtson SH, Gade R, Madsen N, Moeslund TB (2018) Camera calibration for underwater 3d reconstruction based on ray tracing using snell's law. In: 2018 IEEE/CVF Conference on Computer Vision and Pattern Recognition Workshops (CVPRW), pp 1491–14917, DOI 10.1109/CVPRW.2018.00190

Pless R (2003) Using many cameras as one. In: Proceedings of the IEEE Computer Society Conference on Computer Vision and Pattern Recognition, IEEE, vol 2, pp II–587

Rizzini DL, Kallasi F, Oleari F, Caselli S (2015) Investigation of vision-based underwater object detection with multiple datasets. International Journal of Advanced Robotic Systems 12(6):77

Sadowski EM, Schmidt AR, Seyfullah LJ, Solórzano-Kraemer MM, Neumann C, Perrichot V, Hamann C, Milke R, Nascimbene PC (2021) Conservation, preparation and imaging of diverse ambers and their inclusions. Earth-Science Reviews 220:103653

Schonberger JL, Frahm JM (2016) Structure-from-motion revisited. In: Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, pp 4104–4113

Schönberger JL, Frahm JM (2016) Structure-from-motion revisited. In: Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition

Stewenius H, Nistér D, Oskarsson M, Åström K (2005) Solutions to minimal generalized relative pose problems. In: OMNIVIS 2005

Sturm P (2005) Multi-view geometry for general camera models. In: 2005 IEEE Computer Society Conference on Computer Vision and Pattern Recognition (CVPR’05), IEEE, vol 1, pp 206–212

Sweeney C, Fragoso V, Höllerer T, Turk M (2014) gdls: A scalable solution to the generalized pose and scale problem. In: Proceedings of the European Conference on Computer Vision, Springer, pp 16–31

Telem G, Filin S (2010) Photogrammetric modeling of underwater environments. ISPRS journal of photogrammetry and remote sensing 65(5):433–444

Treibitz T, Schechner Y, Kunz C, Singh H (2011) Flat refractive geometry. IEEE transactions on pattern analysis and machine intelligence 34(1):51–65

Ventura J, Arth C, Reitmayr G, Schmalstieg D (2014) A minimal solution to the generalized pose-and-scale problem. In: Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, pp 422–429

Ventura J, Arth C, Lepetit V (2015) An efficient minimal solution for multi-camera motion. In: Proceedings of the IEEE International Conference on Computer Vision, pp 747–755

Xiong J, Heidrich W (2021a) In-the-wild single camera 3d reconstruction through moving water surfaces. In: Proceedings of the IEEE/CVF International Conference on Computer Vision, pp 12558–12567

Xiong J, Heidrich W (2021b) In-the-wild single camera 3d reconstruction through moving water surfaces. In: Proceedings of the IEEE/CVF International Conference on Computer Vision (ICCV), pp 12558–12567

Zhang P, Wu Z, Wang J, Kong S, Tan M, Yu J (2021) An open-source, fiducial-based, underwater stereo visual-inertial localization method with refraction correction. In: 2021 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS), pp 4331–4336, DOI 10.1109/IROS51168.2021.9636198