# Compensating over- and underexposure in optical target pose determination

E. Cledat, M. Rufener, D.A. Cucci

► To cite this version:

E. Cledat, M. Rufener, D.A. Cucci. Compensating over- and underexposure in optical target pose determination. Pattern Recognition, 2021, 116, pp.107930. 10.1016/j.patcog.2021.107930. hal-04424558

HAL Id: hal-04424558

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

Submitted on 29 Jan 2024

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_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Compensating over- and underexposure in optical target pose determination

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

E. Cledat, M. Rufener, D.A. Cucci*

EPFL ENAC TOPO Bâtiment GC - Station 18 1015 Lausanne, Switzerland

## ARTICLE INFO

**Article history:**
Received 9 April 2019
Revised 20 October 2020
Accepted 28 February 2021
Available online 5 March 2021

**Keywords:**
Optical target
Target orientation
Image processing algorithm
Geometry
Ellipse fitting
Computer vision
Overexposure
Exposure compensation
Resection

## ABSTRACT

Optical coded targets allow to determine the relative pose of a camera, on a metric scale, from one image only. Furthermore, they are easily and efficiently detected, opening to a wide range of applications in robotics and computer vision. In this work we describe the effect of pixel saturation and non-ideal lens Point Spread Function, causing the apparent position of the corners and the edges of the target to change as a function of the camera exposure time. This effect, which we call exposure bias, is frequent in over- or underexposed images and introduces a systematic error in the estimated camera pose. We propose an algorithm that is able to estimate and correct for the exposure bias exploiting specific geometric features of a common target design based on concentric circles. Through rigorous laboratory experiments carried out in a highly controlled environment, we demonstrate that the proposed algorithm is seven times more precise and three times more accurate in the target distance estimation than the algorithms available in the literature.

© 2021 Elsevier Ltd. All rights reserved.

## 1. Introduction

Optical coded targets generally consist of a set of high-contrast geometric features (such as lines, squares or circles), in which a code is embedded to exclude false matches and distinguish multiple targets in a scene. Knowing the physical dimensions of the target and the camera calibration allows us to determine from a single image the metric pose of a target relative to the camera. Several target designs and detection algorithms have been proposed in the literature, such as ARToolKit [1], AprilTags [2], and ArUco [3], which have been successfully employed in multiple applications, such as camera calibration, photogrammetry, augmented reality, machine vision and robotics in general.

Target detection algorithms establish correspondences between the geometric features in object space and their projections on the image plane. From these correspondences, the pose of the camera relative to the target can be inferred. For instance, the main target feature in ArUco is a black square on a white background; its four corners are localized by a sub-pixel corner detector and the camera pose is determined by solving the Perspective-n-Points (PnP) problem [4].

The accuracy of the camera pose depends on how precisely the critical target features can be measured on the image. A major source of bias is the exposure time. Indeed, the apparent position of edges and corners may be different in long-exposed images because of the combined effect of a non-ideal lens Point Spread Function (PSF) and pixel saturation. In Fig. 1, we compare two images of an ArUco target taken with two different exposure times: the apparent size of the target is different in the two images and the target on the right seems to be farther away from the camera. We refer to this effect as *exposure bias*, i.e., the systematic error in the estimated camera pose that is a function of the exposure time.

Since saturation in over- or underexposed images leads to a loss of information, it seems difficult to recover the true position of an image feature, such as a corner or an edge, just by using the intensity of pixels in its vicinity. Therefore, image exposure has to be optimized to avoid saturation. However, this is difficult in many cases because of unpredictable light conditions, variable scene illumination as functions of time and space, and reflections on the target surface. Moreover, auto-exposure algorithms are typically implemented at sensor level to optimize either the overall scene exposure or the exposure of a predefined area, which does not guarantee an optimal exposure of the target. As an example, in [5], an UAV follows a terrestrial mapping vehicle that has an optical target mounted on the roof. The pose computed thanks to the optical target have been proved essential to increase the accuracy of reconstructed 3D models by means of photogrammetry. In such outdoor

<footnote>* Corresponding author.
E-mail addresses: emmanuel.cledat@alumni.epfl.ch (E. Cledat), marc.rufener@alumni.epfl.ch (M. Rufener), davide.cucci@epfl.ch (D.A. Cucci).</footnote>
<figure><img src="image_3.png" /><figcaption>**Fig. 1.** Two images of an ArUco target taken from the same camera pose with different exposure times, the underexposed image is on the left and the overexposed one on the right. The brightness of the underexposed image has been enhanced. The difference between the apparent location of the top and bottom edge is highlighted by red and blue lines. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>**Fig. 2.** The proposed target design. Overexposing the target results in a reduction of the apparent thickness of the black ring (purple arrow), while increasing the camera distance leads to a uniform scaling of the whole target. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

conditions, it is difficult to optimise the illumination or exposure of the camera to obtain at the same time an optimal, unsaturated picture of the target, and an optimal illumination in the rest of the scene. This also holds in real-world indoor navigation scenarios, e.g., warehouses, where the light conditions are not uniform and difficult to control. An algorithm to optimize the camera exposure time in real-time was presented by the authors in [6], however this would allow to optimize the target exposure only, possibly compromising the usability of the images for other applications, e.g., in mapping.

In this work, we introduce a new optical target design and its corresponding pose determination algorithm that are not affected by either the exposure time or the illumination. The proposed target design is shown in Fig. 2 and it consists of a black ring over a white background; inside the ring, there are small white dots in which a code is embedded. The same pattern is repeated in a fractal fashion to allow for multi-scale detection (e.g., from very far or very near). In this design, the over- or underexposure affects the apparent thickness of the black ring. However, this effect is orthogonal to the one that we obtain by moving the camera farther or closer to the target, leading to a uniform scaling of the whole target and preserving the proportions of every feature (see the purple and red arrows in Fig. 2). Hence, we can design an algorithm that exploits this effect to estimate and correct the exposure bias in contour positions.

Many studies have investigated how to determine the camera pose from concentric circular features. For instance, in [7], a

method for camera calibration based on concentric circles was presented, and another formulation was proposed in [8]. The geometric and algebraic constraints related to the projection of concentric circles were discussed in [9]. In [10], the authors addressed the issue of how to accurately locate the center of two concentric circles based on their projections (non-concentric ellipses). However, all these works assume that the apparent location of the circles on the image plane is not biased by the exposure. Here, we propose a new algorithm for camera pose determination that exploits the specific geometry of the target to estimate and compensate for any exposure bias in the localization of edges.

This work is organized as follows. In Section 2, we introduce the image processing pipeline for our target design. Each step of the pipeline will be detailed in later sections: in Section 3 and Section 4 we develop a method to determine with sub-pixel precision the apparent image location of the two concentric circles of the target, taking into account a non-ideal lens PSF. In Section 5, we determine the camera pose by using the position of a set of edge points and assuming that those are the projections of two concentric circles of known diameter. In Section 6, we modify such method to compensate for the exposure bias, which we link to the deviations in the apparent thickness of the black ring. In Section 7, we validate our exposure compensation algorithm through acquired ground truth data, proving its illumination invariance and the orthogonality between the distance effect and the exposure bias effect.

## 2. Overview of the image processing pipeline

In this section, we summarize the steps required to determine the position and the orientation of the proposed target relative to the camera.

1. The target is located in the image and its main features, the inner and the outer circle, are coarsely identified. First, a pixel-level edge detection algorithm (e.g., Canny [11]) is applied, then the elliptical contours are clustered in concentric couples. For each couple, the distinctive target code is searched (see again Fig. 2), and the two contours for which the code test succeeds are selected as the outer and inner target circle. This step is described in detail in a previous work of the authors [10].

2. Two ellipses are fitted on the selected contours as in [12]. Then, by using the determined parameters of the ellipses, a set of normal directions $\vec{n}_i$ is analytically computed for a set of equally spaced points $P'_{0,i}$ on the ellipses.

3. The apparent position of the edge along each $\vec{n}_i$ is refined by employing the edge formation model and the sub-pixel edge detector introduced in Section 3 and 4. The result is a set of refined edge points $\vec{P}'_i$.

4. The position and the orientation of the camera relative to the target are computed from $\vec{P}'_i$, without or with exposure compensation, as discussed in Section 5 and 6.

## 3. Edge formation model

In this section, we derive a closed-form model for the light intensity in the vicinity of black/white transitions (i.e., edges) of the target, as measured from the imaging sensor. This model will be used in the next section to construct a rigorous sub-pixel edge detector.

Consider an arbitrary 3D point $P$ on a black/white edge of the target and refer to Fig. 3. Each point on the target maps to a point on the image plane via a function $\Pi(\cdot) : \mathbb{R}^3 \rightarrow \mathbb{R}^2$, so that $P' = \Pi(P)$. $\Pi(\cdot)$ accounts for the camera pose relative to the target, the lens distortion, and the projection onto the image plane. Let now $\vec{n}$ be a unitary vector orthogonal to the edge at $P'$ and
<figure><img src="image_5.png" /><figcaption>Fig. 3. A 3D edge point P maps to its unknown projection P' on the image plane via Π(·). ⃗n is defined as orthogonal to the edge at P', while ⃗m is parallel to it.</figcaption></figure>

directed towards the white side of the edge. We define a refer-
ence frame with axes ⃗n and ⃗m, such that ⃗n ⊥ ⃗m, and with ori-
gin close to P'.

If no significant blur is introduced by the lens system (i.e., the PSF is ideal), the black/white transition occurs sharply on the image plane. More precisely, in the vicinity of P', the light intensity i at coordinates (n, m) on the image plane is independent of m and it is given by

$$
i(m, n) = i(n) = a \mathbb{1}(n - n_{p'}) + b, \quad (1)
$$

where: a and b are two constants such that b is the light intensity
of black areas and (b + a) is the one of white areas, n_{p'} is the n
coordinate of P', and 1(·) is the unit step function.

However, black/white transitions are never sharp since even
properly focused lens systems always introduce a certain amount
of blurring. Thus, the light intensity at each position of the image
plane can be obtained by the convolution of the intensity func-
tion i(m, n) with the lens PSF. We have taken into consideration
only zero-mean Gaussian PSFs as they are widely used to model
lens systems of common cameras. These PSFs are characterized by
a two-dimensional covariance matrix Σ that encodes the spread
of a point light source on the image plane due to a non-ideal lens
system. The light intensity mediated by the lens system is given
by:

$$
\begin{align*}
\mathrm{id}(m, n) &= i(n, m) * \mathrm{psf}(m, n) = i(n, m) * \mathrm{psf}(m, n) \\
&= \int \int_{(\tau_m, \tau_n) \in \mathbb{R}^2} i(n - \tau_n) \mathrm{psf}(\tau_m, \tau_n) \, d\tau_m \, d\tau_n \\
&= \int_{\tau_n \in \mathbb{R}} i(n - \tau_n) \phi(\tau_n; \sigma) \, d\tau_m \, d\tau_n, \tag{2}
\end{align*}
$$

where $\phi(\cdot, \sigma)$ is the one-dimensional Gaussian probability density function with variance $\sigma^2 = \Sigma_{nn}$. Here, we have regarded the light intensity $i(m, n)$ as independent of $m$ in a local neighbourhood of $P'$ because $\vec{m}$ is parallel to the edge, and we have considered the integral in $d\tau_m$ as corresponding to the marginalization of $m$.

By further expanding Eq. (2), we obtain:

$$
\begin{align}
\mathrm{id}(n) = i(n) * \mathrm{pfs}(n) &= \int_{-\infty}^{\infty} (a \mathbb{1}(n - n_{p'} - \tau) + b) \frac{1}{\sigma\sqrt{2\pi}} e^{-\frac{\tau^2}{2\sigma^2}} d\tau \nonumber \\
&= a \int_{-\infty}^{n-n_{p'}} \frac{1}{\sigma\sqrt{2\pi}} e^{-\frac{\tau^2}{2\sigma^2}} d\tau + b = \frac{a}{2} \left[ 1 + \operatorname{erf}\left(\frac{n-n_{p'}}{\sqrt{2}\sigma}\right) \right] + b. \tag{3}
\end{align}
$$

Eq. (3) gives a closed-form expression of the expected light in-
tensity at continuous coordinates on the image plane. It involves
four parameters: i) $n_{p'}$, the actual edge position; ii and iii) a and b,
which encode the light intensity associated with black and white
areas; iv) $\sigma$, the component of the lens PSF in the $\vec{n}$ direction.

In the next section, we show how we can obtain a sub-pixel
estimate of the image coordinates of P' by fitting the model in
Eq. (3) to the pixel intensities in the image.

4. Subpixel edge position

In the following section, we consider the problem of determin-
ing the sub-pixel location of P' given that an initial guess, P'_0, and
the vector orthogonal to the edge and passing trough P'_0, ⃗n, are
known. These are obtained in previous image processing steps. For
example, in [10], the authors determine the pixel-level edge points
of circular target features by using the Canny edge detector. Next,
an ellipse is fitted on those as in [12]; from the ellipse, ⃗n can be
obtained in closed form for each edge point. We also show how
the introduced sub-pixel edge detector is biased if any of the pix-
els in the vicinity of the edge point is saturated, e.g., because of
over- or underexposure.

The sub-pixel edge position P' on ⃗n is obtained by fitting
the edge model in Eq. (3) to the measured pixel intensities along
⃗n via a non-linear least-squares optimization. First, we sample
the image along ⃗n with one pixel spacing, i.e., n_j ∈ [-N, -N +
1, ..., N], where N is the maximum distance from the query point
for which we assume Eq. (3) to hold (e.g., N = 3). Note that the
initial guess for the position of P', P'_0 lies at n = 0. For each n_j,
the image intensity ̂i(n_j) is obtained interpolating the image pix-
els, e.g., by means of bilinear interpolation. Given ̂i(n_i), an estimate
for the location of the edge point along ⃗n, ̂n_{p'}, is calculated by
solving the following least-squares optimization problem:

$$
[\hat{n}_{p'}, \hat{a}, \hat{b}, \hat{\sigma}] = \arg \min_{n_{p'} \text{, } a, b, \sigma} \sum_{j=-N}^{N} (i(n_j) - \tilde{i}(n_j))^2. \quad (4)
$$

This is a non-linear curve-fitting problem that can be solved with
the Levenberg-Marquardt algorithm if the initial guess of every un-
known is sufficiently close. An effective initialization is given by
$n_{p'} = 0$, $a = \tilde{i}(n_N) - \tilde{i}(n_{-N})$, $b = \tilde{i}(n_{-N})$, and $\sigma = 1$. The curve fitting
problem also yields estimates for $a$ and $b$ (the intensity of black
and white areas in the image plane). Note that $\sigma$ is the compo-
nent of the PSF in the direction of $\vec{n}$ and that PSF is stable across
multiple frames and consistent for different areas of the image.

The estimation problem in Eq. (4) is solved for each edge point $P'_{0,i}$, obtaining a set of refined estimates of the edge location along each $\vec{n}_i$. These 1D coordinates are transformed back in the image reference frame and build the set of image observations $\hat{b}'_i$ that we will employ in later sections to estimate the camera pose relative to the target.

We now show that the estimate ˆn_{p'} is biased if some image intensities ˜i(n_j) are saturated. This is common in over- or under-exposed images. Fig. 4 depicts the pixel intensity in an image with an ideal exposure (Fig. 4a) and in an overexposed image (Fig. 4b). In the ideal exposure case, no saturation occurs and the continuous intensity function i(n) (in blue) can be accurately recovered from the samples ˜i(n_j) (in red). In this case, ˆn_{p'} = n_{p'}. On the other hand, in the overexposed image, i(n) is saturated for j > 1. Thus, the estimated intensity function ˆi(n) (in yellow) differs from i(n) in the whole saturated area, resulting in ˆn_{p'} ≠ n_{p'}. Note that the estimated intensity function is close to the true intensity in the non-saturated area (i.e., for n < 2).

Fig. 5 shows a numerical study of the bias in sub-pixel edge position as function of σ and overexposure. The overexposure is parameterized as the ratio between the maximum of the intensity function i(n) and the saturation level. The exposure bias is severe in blurry images (high values of σ) but it is significant also in reasonably sharp images (σ ≈ 1). Similar results hold for the underexposure case.
<figure><img src="image_6.png" /><figcaption>Fig. 4. Pixel intensity of a black/white transition along $\vec{n}$. The intensity $i$ is proportional to the amount of light exposing the sensor. $\tilde{i}$ is the measured pixel intensity (subject to overexposure) and $\tilde{i}$ is the estimated from $\tilde{i}$. The estimated edge position $\hat{n}_p$ differs from the true value $n_p$ by a few hundredths of a pixel if the image exposure is ideal. However, in overexposed images the difference can be of multiple pixels, see Fig. 5.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Fig. 5. Exposure bias as a function of the overexposure ratio for different values of $\sigma$.</figcaption></figure>

To give a measure of the error in the camera pose determination due to the exposure bias, we introduce an error $\Delta d$ when we determine the target distance $d$ by comparing its known size $2r$ with the size measured on the image plane $2r'$ (which is subject to the exposure bias $\gamma$). When a target is parallel to the image plane, this error is given by $\Delta d = -\frac{d\gamma}{r' + \gamma} \approx -d\gamma$. $\Delta d$ is proportional to the distance from the target, and it is approximately proportional to the relative exposure bias (i.e., exposure bias normalized by the target radius). For example, a 5 cm target imaged from 3 m, with focal length $f = 2900$ px, has a diameter of 48 px in the image. In this case, an exposure bias of one pixel leads to an error of 13 cm ($\approx 5\%$) in the distance estimation. Such error affects any target detection library that relies on the image positions of certain target features: if these are biased by exposure, the resulting distance estimate will be biased as well. To show this, we compute the distance to a 5 cm square target (same size as our inner circular target) with the ArUcO library, from a distance of approximately 3 m. The results are shown in Fig. 6 and are consistent with the aforementioned theoretical derivations.

In the next sections, we will explain how specific features of our target design (Fig. 2) allow estimating and correcting for the exposure bias.

## 5. Camera pose determination

We present a new rigorous approach to determine the camera pose starting from a set of edge points being the projection of the two concentric circles composing the target (see Fig. 2). In Section 6, we will extend this approach to estimate and compensate for the exposure bias.

After applying the sub-pixel detection algorithm introduced in Section 4, a set of edge points $\hat{P}'_i$ is available. These are the projections of object points lying on two concentric circles of known radius. The camera pose is determined by minimizing the projection error associated with each of these points.

The edge points are given in homogeneous coordinates $\hat{P}'_i = [x_{\hat{P}'_i}, y_{\hat{P}'_i}, 1]^T$ and they will be used as input for the algorithm presented in this section. First, $\hat{P}'_i$ is transformed into a unit-less vector $\tilde{P}'_i$ by employing the camera calibration matrix:

$$ \tilde{P}'_i = \begin{pmatrix} f & 0 & c_x \\ 0 & f & c_y \\ 0 & 0 & 1 \end{pmatrix}^{-1} \hat{P}'_i \quad (5) $$

where $f$ is the principal distance of the camera, and $c_x$ and $c_y$ are the $x$ and $y$ coordinate of the principal point, respectively. Eq. (5) can also be modified to account for lens distortion, as described in [13]. Note that the intrinsic camera calibration must be known a priori because the focal length $f$ is correlated with the depth of the target, and they cannot be simultaneously determined with the approach presented in this work.

Let the target reference frame be placed at the center of the two concentric circles of the target with the Z axis normal to the target plane and pointing away from it. Let $T$ be the position of the optical center of the camera in this frame, and $R$ the rotation from the camera to this frame (Fig. 7). Since $\hat{P}'_i$ can be the projection of a point belonging either to the outer or to the inner circle of the target, we define the function $r(i) : \mathbb{N} \to \mathbb{R}_+$ that provides the radius of the circle to which $\hat{P}'_i$ belongs. It is worth noting that this function allows us to generalize the approach presented below to an arbitrary number of concentric circles.

We know from perspective geometry that, for any object point $P$ and its projection onto the image plane $\tilde{P}'$, given in unit-less homogeneous coordinates, it holds that

$$ \exists \lambda \in \mathbb{R} \mid T + \lambda R \tilde{P}' = P. \quad (6) $$
<figure><img src="image_8.png" /><figcaption>Fig. 6. Exposure bias and distance for the target in Fig. 1 as computed by the ArUco library, as a function of the image exposure time. The exposure bias -in this context- is defined as the (mean) displacement of the Aruco target edges toward black side, in pixels, with respect to the reference position taken for an exposure time of 50 000 µs.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Fig. 7. Geometric elements involved in the pose determination problem. The image plane is depicted in red. The pose of the camera reference frame in relation to the target frame is given by T and R. An object point $P_i$ lying on one of the two concentric circles of the target with angle $\theta_i$ projects to $\tilde{P}'_i$ on the image plane. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

$\lambda$ is generally obtained from the third component of Eq. (6) and eliminated.

For any $\tilde{P}'_i$, we know that the corresponding object point lies on a circle of radius $r(i)$ on the plane $Z = 0$. Therefore we have that:

$$ \exists \lambda_i \in \mathbb{R}, \theta_i \in [0, 2\pi) \mid T + \lambda R \tilde{P}'_i = \begin{bmatrix} r(i) \cos \theta_i \\ r(i) \sin \theta_i \\ 0 \end{bmatrix}. \quad (7) $$

$\theta_i$ encodes, for each edge points, its corresponding position on the target circle in object space: it is unknown at this stage since it is not possible, by using only local features of the target, to establish

an unambiguous 2D-to-3D correspondence between object points on the target circles and their projections on the image plane. We eliminate $\lambda$ and $\theta_i$ from Eq. (7), as discussed in detail in Appendix, obtaining a 1D condition for each point $\tilde{P}'_i$.

We can formulate Eq. (7) for each edge point $\tilde{P}'_i$. This gives us a system of equations that can be solved for the desired R and T (e.g., in least-squares sense), provided that enough edge points are available on both circles. However, since we have not established an unambiguous 2D-to-3D correspondence, such system is under-constrained. Indeed, Eq. (7) is still satisfied if we multiply the left side of the equality by an arbitrary rotation matrix $\Omega_Z$ around the Z axis. This issue can be solved by constraining the optical center of the camera to the Y = 0 plane. This is achieved as follows: let

$$ \begin{cases} \tilde{T} = \Omega_Z T \\ \tilde{R} = \Omega_Z R' \end{cases} \quad (8) $$

such that $\tilde{T} = [\tilde{X}, 0, \tilde{Z}]$. This requires that

$$ \Omega_Z = \begin{bmatrix} \frac{T_X}{\sqrt{T_X^2 + T_Y^2}} & \frac{T_Y}{\sqrt{T_X^2 + T_Y^2}} & 0 \\ -\frac{T_Y}{\sqrt{T_X^2 + T_Y^2}} & \frac{T_X}{\sqrt{T_X^2 + T_Y^2}} & 0 \\ 0 & 0 & 1 \end{bmatrix}. \quad (9) $$

We can now replace T and R with $\tilde{T}$ and $\tilde{R}$ in Eq. (7) and rearrange it in the form

$$ f_i(\tilde{T}, \tilde{R}, \tilde{P}'_i) = 0. \quad (10) $$

This gives us one condition $f_i$ per observation $\tilde{P}'_i$, providing a system of equations that can be solved for $\tilde{T}$ and $\tilde{R}$, e.g., in the least-squares sense, with the Gauss-Newton or Levenberg-Marquardt algorithm. As $\tilde{R}$ belongs to the special orthonormal group SO(3), we rely on the manifold encapsulation technique described in [14] to deal with non-Euclidean unknowns in least-squares estimation.

As already mentioned, the target can rotate around the Z axis without changing the geometry of the problem (refer to Fig. 7). This is because we cannot establish unique 3D-to-2D correspondences between object points on the target circles and their image projections. Indeed, $\theta_i$ is not known in Eq. (7). To solve this issue,
<figure><img src="image_10.png" /><figcaption>Simple concentric circle target</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>Pixel value on diameter line</figcaption></figure>

<caption>**Fig. 8.** Pixel intensity along the diameter of the target for an underexposed (top) and an overexposed (bottom) image. The brightness of the top image has been enhanced. In the underexposed image, the black/white transition is shifted by a few pixels towards the black area, which is highlighted by red and blue lines in the plots on the right. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)</caption>

we constrain the camera to the $Y = 0$ plane. However, this is not a limiting problem since the position of the target relative to the camera, $-\hat{R}^T \hat{T}$, can still be determined.

In this section, we have described a method to determine the pose of the camera relative to the target by using a set of edge points being the projection of two concentric circles with known radius. In the next section, we will extend this method to take into account the exposure bias, i.e., the condition in which the location of the edge points on the image plane is affected by over- or underexposure.

## 6. Correcting for the exposure bias

In the previous section, we have seen how to determine the position and the orientation of the camera relative to the proposed

target from a set of sub-pixel edge locations obtained as described in Section 4. However, we also saw that these measurements can be biased if the exposure time is not ideal (see Fig. 4). Indeed, in case of over- or underexposure, the apparent position of $\hat{P}'_i$ can be displaced along the line normal to the edge $\tilde{n}_i$. When the image is overexposed, such displacement is towards black areas and vice versa for underexposed images. This effect is shown in Fig. 8, which compares the pixel intensities of two pictures taken with different exposure times from the same camera position and orientation.

Recall that, as defined in Section 3, $\tilde{n}_i$ is a unitary vector on the image plane, perpendicular to the edge, and directed towards the white side of the edge. Let $\gamma_i \in \mathbb{R}$ be the shift of $\tilde{P}'_i$ along the vector $\tilde{n}_i$. If the target has a uniform illumination, we can assume $\gamma_i = \gamma$, $\forall i$. The new parameter $\gamma$ can be estimated along the geo-
<figure><img src="image_12.png" /><figcaption>Fig. 9. Reprojection error for the considered edge points $\hat{P}'_i$ in an image without (on the left) and with (on the right) exposure compensation.</figcaption></figure>

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

metric parameters $\hat{T}$ and $\hat{R}$. In order to do this, we substitute $\hat{P}'_i$ in Eq. (7) with an observation corrected for the exposure: $\hat{P}'_i + \gamma \hat{n}_i$. Hence:

$$ \exists \lambda_i \in \mathbb{R}, \theta_i \in [0, 2\pi) \mid T + \lambda R(\hat{P}'_i + \gamma \hat{n}_i) = \begin{bmatrix} r(i) \cos \theta_i \\ r(i) \sin \theta_i \\ 0 \end{bmatrix}. \quad (11) $$

The three conditions of this system can be simplified to one condition in the form $f_i(\hat{T}, \hat{R}, \gamma, \hat{P}'_i) = 0$ for each edge measurement $\hat{P}'_i$ by eliminating $\lambda$ and $\theta_i$, as discussed in the Appendix.

Fig. 9 shows the residuals associated with each $f_i$ in an overexposed image with (on the right) and without (on the left) correction for exposure bias. In the case of overexposure, the thickness of the black ring is smaller then expected, which translates in systematic residuals that cannot be compensated by any assignment of $\hat{R}$ and $\hat{T}$ (the residuals on the inner circle points are bigger as those points are outnumbered by the ones on the outer circle). On the other hand, such systematic effect in the residuals is completely captured once the estimation of the $\gamma$ parameter is enabled. The impact of this on the accuracy of the camera pose determination with respect to the target will be investigated in the next section.

<figure><img src="image_14.png" /><figcaption>Fig. 10. Schematic representation of the set-up to obtain the reference positions of the target and the camera. The theodolites network in B, B' and B'' measures the position of the camera (A), the position of several ArUco targets (C), and the position and orientation of the proposed target (D). The camera A is oriented from the theodolites and from image observation of the ArUco targets, C.</figcaption></figure>

## 7. Experimental evaluation

The goal of our method is to determine the position of a camera relative to a target compensating for any possible bias caused by over- or underexposure. We validated it in a highly controlled environment: refer to Fig. 10. Both the camera A and the target D were fixed on a rigid tripod, and a set of photos with different exposure times was taken to simulate normal or difficult light conditions. The same set-up was repeated for 20 different positions of the target, and, for each of them, the relative position of the camera was determined i) by means of theodolites to compute a ground truth positioning, and ii) with several Computer-Vision algorithms, including the ones presented in this paper. The results are then compared focusing on the estimated distance from the camera to the target, since the effects of over- and underexposure mainly translates in a bias in the apparent target distance.

The concentric circle target D was printed on a 50 × 50 cm white plastic rigid board (the target design is shown in Fig. 2), and the theodolites measured four crosses printed on the four corners of the board (Fig. 10).

This experimental set-up allows computing the position of the camera A relative to the target D. By using the images taken by the

camera A, we validated the method presented in Sections 5 and 6 and determined the camera position. Since the pose determination depends on the exposure time, a set of 117 photos was taken with different exposure times ranging from 10 ms to 300 ms (such high exposure times depends on the rather dark indoor setup used for the experiment). We applied four computer-vision algorithms to determine the position of A relative to D, considering both the outer target (diameter: 25cm) and the inner target (diameter: 5cm). We first applied the PnP algorithm [15] by using the centroids of the five white dots composing the target code (Fig. 2); second, we used an algorithm based on separated ellipse fitting [10]; third, we run the joined ellipse adjustment described in Section 5; lastly, we corrected the joined ellipse adjustment for the exposure bias as described in Section 6.

We note that it is not possible to apply a state-of-the-art target detection algorithm, such as ArUco, to our target design because of the radical differences in the target shape and figures. However, in Section 4 we showed theoretically and experimentally that ArUco, and any other library that is not based on exposure invariant image features, would be affected by exposure bias. State-of-the-art
<figure><img src="image_15.png" /><figcaption>**Fig. 11.** Target distance as a function of exposure time obtained with different algorithms with respect to the reference. It is possible to see that the distance obtained with exposure compensation, as in **Section 6**, is practically independent with respect to exposure time.</figcaption></figure>

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

target detection algorithms solve a PnP problem to recover the target pose from the measured image features, e.g., the four corners of a squared target in ArUcO, see [3], **Section 4**. We include the PnP algorithm in the experimental evaluation to show that if exposure-invariant image features are employed, such as the centroids of the five white dots composing our target code (Figure 2), the exposure bias disappears, however leading to a more noisy distance estimate with respect to our solution, as we will discuss in details in the following.

Fig. 11 shows the computed target distance with the four different algorithms for a given camera position and for each exposure time. First, when comparing the obtained results with the reference distance, we notice that the errors associated with the inner target are about five times bigger than the errors for the outer target. This is expected, as the inner target is five times smaller than the outer one. Second, we notice that when the image is not over-exposed (e.g., exposure time less than 50 ms, no algorithm exhibits a bias which can be correlated with the exposure time. However, if the image was overexposed (e.g., exposure time greater than 50 ms), all the algorithms without exposure compensation, with the exception of PnP, i.e., ([10] and **Section 5**) show a deviation from the reference distance that increased with the exposure time. On the other hand, when the exposure bias $\gamma$ was computed together with the other parameters (**Section 6**), the results became independent from the exposure. Regarding the PnP algorithm, in this particular case it is not affected by exposure time as the position of the white dots, used as an input for this algorithm, does not change as a function of the exposure. This would not be the case, for example, in the case corner points are used, as in ArUcO, as shown in Figure 6. However, since only five points are available, the distance obtained from PnP is more noisy if compared with the other algorithm employing hundreds of points on the two circles.

The same experiment was repeated for different positions of the target D. The position of A, B, B', B'' and C remained unchanged while we modified the position of D 20 times. Two statistical indicators were calculated to assess the quality of the results: the Standard deviation and the RMS of true error (Fig. 12). The Standard deviation is determined by comparing the results of each exposure time to their mean $\bar{d}$; it describes the dispersion of the computed values but it does not assess the influence of a potential bias.

<caption>**Table 1**</caption>

Geometric mean of the relative error for every position. OT stands for Outer target and IT stands for Inner target.

<table><thead><tr><th rowspan="3"></th><th colspan="2">Geometric mean of:</th></tr><tr><th>relative<br>standard deviation</th><th>relative<br>RMS of true error</th></tr><tr><th><sup>1</sup>&frasl;<sub>2</sub> &sum;<sub>i=1</sub><sup>n</sup> (d<sub>i</sub> - d&#773;)<sup>2</sup><br><u>d</u><sup>d</sup></th><th><sup>1</sup>&frasl;<sub>2</sub> &sum;<sub>i=1</sub><sup>n</sup> (d<sub>i</sub> - d&#773;)<sup>2</sup><br><u>d</u><sup>d</sup></th></tr></thead><tbody><tr><td>OT Without exposure compensation</td><td>0.49 %</td><td>1.01 %</td></tr><tr><td>With exposure compensation</td><td>0.05 %</td><td>0.37 %</td></tr><tr><td>Cucci, 2016 [8]</td><td>0.66 %</td><td>1.17 %</td></tr><tr><td>PnP</td><td>0.10 %</td><td>0.15 %</td></tr><tr><td>IT Without exposure compensation</td><td>2.79 %</td><td>4.30 %</td></tr><tr><td>With exposure compensation</td><td>0.44 %</td><td>1.16 %</td></tr><tr><td>Cucci, 2016 [8]</td><td>4.16 %</td><td>5.73 %</td></tr><tr><td>PnP</td><td>1.76 %</td><td>1.86 %</td></tr></tbody></table>

To address this issue, we also computed the RMS of true error based on the comparison between the obtained results and their ground truth value $\bar{d}$.

The summary statistics are shown in **Table 1**. The errors in Fig. 12 are normalized by the distance to the target and aggregated by geometric mean. These results show that the exposure compensation is more than seven times more precise, and more than three times more accurate than ellipse fitting without exposure compensation. However, since the position of the white dots is not affected by the exposure, the distance calculated by the proposed method is only about two times more precise than the one assessed by the PnP, and it is approximately as accurate.

We validated the method described in **Section 6** also by assessing the independence between the geometric parameters and the exposure bias parameter $\gamma$ by calculating the correlation factor between the computed distance $d$ and $\gamma$. The covariance matrix of $\bar{X}$, $\bar{Z}$, and $\gamma$ is computed after solving the least-squares estimation problem formulated in **Section 6**. The covariance matrix of $d$ and $\gamma$ is determined by variance propagation as follows:

$$ \begin{bmatrix} \sigma_d^2 & \sigma_{d\gamma} \\ \sigma_{d\gamma} & \sigma_{\gamma}^2 \end{bmatrix} = F \begin{bmatrix} \sigma_{\bar{X}}^2 & \sigma_{\bar{X}\bar{Z}} & \sigma_{\bar{X}\gamma} \\ \sigma_{\bar{X}\bar{Z}} & \sigma_{\bar{Z}}^2 & \sigma_{\bar{Z}\gamma} \\ \sigma_{\bar{X}\gamma} & \sigma_{\bar{Z}\gamma} & \sigma_{\gamma}^2 \end{bmatrix} F^T, \quad (12) $$
<figure><img src="image_17.png" /><figcaption>Fig. 12. Distance error as a function of the true target distance for different pose determination algorithms.</figcaption></figure>

where the *F* matrix is given below:

$$F = \begin{bmatrix} \frac{\tilde{x}}{\sqrt{\tilde{x}^2 + \tilde{z}^2}} & \frac{\tilde{z}}{\sqrt{\tilde{x}^2 + \tilde{z}^2}} & 0 \\ 0 & 0 & 1 \end{bmatrix}. \quad (13)$$

The correlation between *d* and *γ* was computed by normalizing the non-diagonal term of the covariance matrix between *d* and *γ* (given by 12) by the square root of their variances: $\sigma_{d\gamma}/(\sigma_d \cdot \sigma_\gamma)$. The typical value was around 5%, proving the orthogonality between the determination of the distance from the camera to the target and the determination of the exposure effect on the position of the ellipse edge (as suggested by Fig. 2).

## 8. Conclusions

We have introduced a new method to rigorously fit the image of a concentric circle target to determine the pose of the camera that sees it. Moreover, we have addressed the problem of the exposure that displaces black and white edges, moving them towards the white direction and leading to systematic effects. Then, we have proposed a way to estimate this displacement and reduce the error.

We performed experimental evaluations of the algorithm in controlled laboratory conditions: We measured the reference target position with theodolites with millimeter-accuracy and captured several images with different exposure times; we then compared the measured distance with the distances calculated by our method at every exposure time. Our method led to highly accurate estimates of the target position with a mean error of 0.37% of the

distance on the 25 cm diameter target, the accuracy of the measurements was also not affected by changes in the illumination. We showed the illumination invariance property of the novel algorithm for a single target position by plotting the target distance estimated with the algorithm in [10] and the illumination invariant one. Finally, we observed that the *γ* parameter is basically uncorrelated (< 5%) with the target position parameters.

## Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

## Acknowledgments

This research has been partially funded by the European Community Horizon 2020 Programme under Grant No. 641518: "EGNOS-GPS/GALILEO-based high-resolution terrestrial-aerial sensing system" managed by the European GNSS Agency (GSA).

## Appendix A. Derivatives of *f* with respect to parameters and observations

The goal of this appendix is to simplify Eq. (7) and Eq. (11) in order to clarify the least-square optimization problem. For convenience, we set a new vector $\tilde{P}_i$ that is equal to $\tilde{R} \tilde{P}'_i$ when there is no exposure compensation (as in Section 5), or to $\tilde{R}(\tilde{P}'_i + \gamma \tilde{n}_i)$
<figure><img src="image_18.png" /><figcaption>Fig. A1. Calculation flow. At the top, we see the input, the intermediates steps, and the output; at the bottom, the formulae used to calculate the parameters.</figcaption></figure>

when there is exposure compensation (as in Section 6). Thus, both Eq. (7) and Eq. (11) can be simplified as below (where the vector $\tilde{P}_{XY_i}$ is composed of the two first components of $\tilde{P}_i$, and the vector $\tilde{P}_{Z_i}$ of its last component):

$$ \left\{ \begin{array}{l} \left\| \begin{pmatrix} \tilde{X} \\ 0 \end{pmatrix} + \lambda_i \tilde{P}_{XY_i} \right\| = r(i) \\ \tilde{Z} + \lambda_i \tilde{P}_{Z_i} = 0 \end{array} \right. \qquad (14) $$

This system of two equations can be solved by removing $\lambda_i$. Therefore, for each measurement point, the condition described in Eq. (15) is imposed:

$$ \left\| \begin{pmatrix} \tilde{X} \\ 0 \end{pmatrix} - \frac{\tilde{Z}}{\tilde{P}_{Z_i}} \tilde{P}_{XY_i} \right\|^2 - r(i)^2 = 0. \qquad (15) $$

To solve the non-linear least-squares optimization problem, we need to compute the derivatives of the conditions with respect to the parameter vectors and the observation vectors.

For convenience, we introduce the notations $\tilde{P}$, $\Psi$ and $\mathfrak{N}$:

$$ \tilde{P}_i = \frac{\tilde{P}_{XY_i}}{\tilde{P}_{Z_i}}, \qquad (16) $$

$$ \Psi_i = \begin{pmatrix} \tilde{X} \\ 0 \end{pmatrix} - \tilde{Z} \tilde{P}_i, \qquad (17) $$

$$ \mathfrak{N}_i = \|\Psi_i\|^2. \qquad (18) $$

Hence, Eq. (15) can be re-written as:

$$ \left\{ \begin{matrix} \begin{pmatrix} \tilde{X} \\ 0 \end{pmatrix} - \tilde{Z} \frac{\tilde{P}_{XY_i}}{\tilde{P}_{Z_i}} \tilde{P}_i \\ \Psi_i \end{matrix} \right\}^2 - r(i)^2 = 0. \qquad (19) $$

Fig. A13 sums up the calculation flow to compute a condition such as 15.

Fig. A13 helps to compute the derivatives of 15 with respect to the parameters $\frac{\partial \mathfrak{N}}{\partial \tilde{X}}$, $\frac{\partial \mathfrak{N}}{\partial \tilde{Z}}$, $\frac{\partial \mathfrak{N}}{\partial \tilde{R}}$ or to the observations $\frac{\partial \mathfrak{N}}{\partial \tilde{P}'}$ (here, we have intentionally misused a derivative with respect to a rotation matrix. The rigorous approach is described below).

First, we introduce the computation of the Jacobian matrix of the condition with respect to the observations:

$$ \frac{\partial \mathfrak{N}}{\partial \tilde{P}'} = \frac{\partial \mathfrak{N}}{\partial \Psi} \frac{\partial \Psi}{\partial \tilde{P}} \frac{\partial \tilde{P}}{\partial \tilde{P}'} \frac{\partial \tilde{P}}{\partial \tilde{P}'} \qquad (20) $$

The four (elementary) matrices to compute are derived below:

$$ \frac{\partial \mathfrak{N}}{\partial \Psi} = 2(\Psi_X \quad \Psi_Y), \qquad (21) $$

$$ \frac{\partial \Psi}{\partial \tilde{P}} = -\tilde{Z} I_2, \qquad (22) $$

$$ \frac{\partial \bar{P}}{\partial \bar{P}} = \frac{1}{\bar{P}_Z} \begin{pmatrix} 1 & 0 & -\bar{P}_X/\bar{P}_Z \\ 0 & 1 & -\bar{P}_Y/\bar{P}_Z \end{pmatrix}, \qquad (23) $$

$$ \frac{\partial \bar{P}}{\partial \bar{P}'} = \tilde{R}. \qquad (24) $$

Similarly to the previous calculation, it is possible to compute the derivative with respect to $\tilde{X}$:

$$ \frac{\partial \mathfrak{N}}{\partial \tilde{X}} = \frac{\partial \mathfrak{N}}{\partial \Psi} \frac{\partial \Psi}{\partial \tilde{X}}. \qquad (25) $$

A new matrix needs to be computed:

$$ \frac{\partial \Psi}{\partial \tilde{X}} = \begin{pmatrix} 1 \\ 0 \end{pmatrix}. \qquad (26) $$

Similarly to the previous calculation, it is possible to compute the derivative with respect to $\tilde{Z}$:

$$ \frac{\partial \mathfrak{N}}{\partial \tilde{Z}} = \frac{\partial \mathfrak{N}}{\partial \Psi} \frac{\partial \Psi}{\partial \tilde{Z}}, \qquad (27) $$

and a new matrix needs to be calculated:

$$ \frac{\partial \Psi}{\partial \tilde{Z}} = -\tilde{P}. \qquad (28) $$

The Jacobian matrix with respect to $R$ must use methods to properly handle the rotations, such as the theory of Lie-groups. Here we give a short explanation of how this theory can be applied to the adjustment.

The position of the camera is described by two parameters: $\tilde{X}$ and $\tilde{Z}$. These parameters belong to an Euclidean space, and, for this reason, we can use the standard actualization step:

$$ \underbrace{\begin{pmatrix} \tilde{X} \\ \tilde{Z} \end{pmatrix}}_{\text{Adjusted}} = \underbrace{\begin{pmatrix} \tilde{X} \\ \tilde{Z} \end{pmatrix}}_{\text{Approximated}} + \begin{pmatrix} \delta \tilde{X} \\ \delta \tilde{Z} \end{pmatrix}. \qquad (29) $$

However, the rotations do not belong to an Euclidean space but to the SO(3) rotation group. Therefore, we need to change their actualization step:

$$ \underbrace{\tilde{R}}_{\text{Adjusted}} = \exp([\omega]_{\times}) \underbrace{\tilde{R}}_{\text{Approximated}}, \qquad (30) $$

where $\omega$ does not belong to the non-Euclidean rotation group but to the Euclidean tangent plane of $R$ (see [14]). The derivatives are computed with respect to $\omega$. The exponential function allows converting a value from the tangent plane into the Lie-Group. $[\omega]_\times$ is the Skew-symmetric matrix built from $\omega$. We can now compute the Jacobian matrix of the condition with respect to the rotation:

$$ \frac{\partial \mathfrak{N}}{\partial \omega} = \frac{\partial \mathfrak{N}}{\partial \Psi} \frac{\partial \Psi}{\partial \tilde{P}} \frac{\partial \tilde{P}}{\partial \tilde{P}} \frac{\partial \tilde{P}}{\partial \omega}. \qquad (31) $$

For any 3D vector $\mathbf{x}$, if $\omega$ and $R$ are linked by the Eq. (30), the following formula can be introduced [14]:

$$ \frac{\partial (Rx)}{\partial \omega} = -[Rx]_\times, \qquad (32) $$

which can be applied to the Eq. (31), giving:

$$ \frac{\partial \bar{P}}{\partial \omega} = -[\bar{P}]_\times. \qquad (33) $$

The Jacobian matrices computed above allow running the adjustment with the five geometric parameters of our initial question
(two parameters for the camera position, and three for its rota-
tion). If we want to take into account the exposure compensation
(i.e., if the parameter $\gamma$ is not fixed to 0), we also need to compute
the Jacobian matrix with respect to this parameter:

$$
\frac{\partial \mathbf{\kappa}}{\partial \gamma} = \frac{\partial \mathbf{\kappa}}{\partial \psi} \frac{\partial \psi}{\partial \bar{\mathbf{P}}} \frac{\partial \bar{\mathbf{P}}}{\partial \bar{\mathbf{P}}} \frac{\partial \bar{\mathbf{P}}}{\partial \gamma}. \quad (34)
$$

A new matrix needs to be computed:

$$
\frac{\partial \bar{P}}{\partial \gamma} = \tilde{R} \tilde{n}. \qquad (35)
$$

The Jacobian calculations described in this appendix must be
determined for each condition, i.e., for each measurement. Then, it
is possible to build the matrices **A** and **B**, which are the Jacobian
matrices of the function **f** with respect to parameters and observa-
tions, respectively (where **f** is the concatenation of all conditions
$f_i$ defined in Figure A1). For each iteration of the Gauss-Newton or
the Levenberg-Marquardt algorithm, we can formulate and solve
the following linearised equation:

$$
\left[
\begin{array}{ccccccc}
\frac{\partial \mathbf{\kappa}_1}{\partial \tilde{\mathbf{X}}} & \cdots & \frac{\partial \mathbf{\kappa}_1}{\partial \tilde{\mathbf{Z}}} & \cdots & \frac{\partial \mathbf{\kappa}_1}{\partial \omega_X} & \cdots & \frac{\partial \mathbf{\kappa}_1}{\partial \omega_Z} & \cdots \\
\vdots & \ddots & \vdots & \ddots & \vdots & \ddots & \vdots & \ddots \\
\frac{\partial \mathbf{\kappa}_1}{\partial \tilde{\mathbf{Y}}} & \cdots & \frac{\partial \mathbf{\kappa}_1}{\partial \omega_X} & \cdots & \frac{\partial \mathbf{\kappa}_1}{\partial \omega_Z} & \cdots & \frac{\partial \mathbf{\kappa}_1}{\partial \gamma} & \cdots \\
\hline
\multicolumn{8}{c}{\left\{ \begin{array}{c} \mathbf{\kappa} \\ \mathbf{\kappa} \end{array} \right\}}
\end{array}
\right]
+ \left[
\begin{array}{c}
0 \\
\frac{\partial \mathbf{\kappa}_1}{\partial \tilde{P}'_i} \\
\vdots \\
0
\end{array}
\right]
\underbrace{\left[
\begin{array}{c}
\mathbf{v}_i \\
\vdots \\
\mathbf{w}_i
\end{array}
\right]}_{\mathbf{v}} + \underbrace{\left[
\begin{array}{c}
\mathbf{w}_i \\
\vdots \\
0
\end{array}
\right]}_{\mathbf{w}} = \underbrace{\left[
\begin{array}{c}
\mathbf{\kappa} \\
\mathbf{\kappa} \\
\mathbf{\kappa} \\
\vdots \\
\mathbf{\kappa}
\end{array}
\right]}_{\mathbf{\kappa}}. \quad (36)
$$

The misclosures $w_i$ are computed from the observations $\tilde{p}'_i$
thanks to the condition $f_i$:

$$
w_i = f_i(\tilde{T}, \tilde{R}, \tilde{P}'_i, \gamma). \tag{37}
$$

This problem is solved in a least-square sense (the L2 norm
of the residuals vector $\boldsymbol{v}$ is minimized) by using minimization al-
gorithms, such as the Gauss-Newton or the Levenberg-Marquardt
method. At each step of the algorithm, the actualisation step is
given by Eq. (29) for the position, and by Eq. (30) for the orien-
tation. Moreover, data-snooping is performed to remove conditions
that have high normalized misclosure.

References

[1] D. Wagner, G. Reitmayr, A. Mulloni, T. Drummond, D. Schmalstieg, Pose track-
ing from natural features on mobile phones, in: Proceedings of the 7th
IEEE/ACM International Symposium on Mixed and Augmented Reality, IEEE
Computer Society, 2008, pp. 125–134.

[2] E. Olson, AprilTag: A robust and flexible visual fiducial system, in: Proceedings of the IEEE International Conference on Robotics and Automation (ICRA), IEEE, 2011, pp. 3400–3407.

[3] S. Garrido-Jurado, R. Muoz-Salinas, F. Madrid-Cuevas, M. Marn-Jimnez, Automatic generation and detection of highly reliable fiducial markers under occlusion, Pattern Recognit 47 (6) (2014) 2280–2292, doi:10.1016/j.patcog.2014.01.005.

[4] Y. Wu, Z. Hu, Pnp problem revisited, J Math Imaging Vis 24 (1) (2006) 131–141.

[5] P. Molina, M. Blázquez, D.A. Cucci, I. Colomina, First results of a tandem terres-
trial-unmanned aerial mapkite system with kinematic ground control points
for corridor mapping, Remote Sens (Basel) 9 (1) (2017) 60.

[6] L. Jospin, A. Stoven-Dubois, D.A. Cucci, Photometric long-range positioning of led targets for cooperative navigation in uavs, Drones 3 (3) (2019) 69.

[7] G. Jiang, L. Quan, Detection of concentric circles for camera calibration, in: Computer Vision, 2005. ICCV 2005. Tenth IEEE International Conference on, 1, IEEE, 2005, pp. 333–340.

[8] F. Abad, E. Camahort, R. Vivó, Camera calibration using two concentric circles, in: International Conference Image Analysis and Recognition, Springer, 2004, pp. 688–696.

[9] J.-S. Kim, P. Gurdjos, I.-S. Kweon, Geometric and algebraic constraints of pro-
jected concentric circles and their applications to camera calibration, IEEE
Trans Pattern Anal Mach Intell 27 (4) (2005) 637–642, doi:10.1109/TPAMI.2005.
80.

[10] D.A. Cucci, Accurate optical target pose determination for applications in aerial photogrammetry, ISPRS Annals of Photogrammetry, Remote Sensing and Spatial Information Sciences III-3 (2016) 257–262, doi:10.5194/ isprs-annals-III-3-257-2016.

[11] J. Canny, A computational approach to edge detection, IEEE Trans Pattern Anal
Mach Intell (6) (1986) 679–698.

[12] A.W. Fitzgibbon, R.B. Fisher, et al., A buyer's guide to conic fitting, DAI Research paper (1996).

[13] D.C. Brown, Close-range camera calibration, PHOTOGRAMMETRIC ENGINEERING 37 (8) (1971) 855–866.

[14] H. Strasdat, Local Accuracy and Global Consistency for Efficient Visual SLAM, Imperial College London - Department of Computing, 2012 Ph.D. thesis.

[15] V. Lepetit, F. Moreno-Noguer, P. Fua, Epnp: an accurate o(n) solution to the pnp problem, Int J Comput Vis 81 (2) (2008) 155, doi:10.1007/s11263-008-0152-6.

Emmanuel Cledat received his master degree in geomatic engineering from Institut national des sciences appliquées (INSA), Strasbourg, France in 2015. His master thesis focused on 3D modeling of complex industrial environments using terrestrial laser scanning and was run in collaboration with Électricité de France (EDF). He received his Ph.D. from École polytechnique fédérale de Lausanne (EPFL), Switzerland, for his thesis under the supervision of Prof. Jan Skaloud. He is currently associate professor at Institut géographique national (IGN) and École Nationale des Sciences Géographiques (ENSG). His current research interests are image processing, computer vision, photogrammetry, laser scanning and geographic information systems.

Marc Rufener received his bachelor and master degree in environmental engineer-
ing from École polytechnique fédérale de Lausanne (EPFL), Switzerland, in 2015 and
2017. His master thesis involved the optimization of a source-receptor air quality
model to simulate the effects of emission reductions. He is currently database man-
ager at the canton administration of Neuchâtel, Switzerland. His research interests
span spatial data analysis, 3D modeling and geomatics.

Davide A. Cucci received his bachelor and master degrees in computer engineer-
ing from Politecnico di Milano, Italy, in 2008 and 2011. From the same university,
he obtained the Ph.D. degree on 2014. His thesis covered the development of multi-
sensor fusion methods for the real-time navigation and mapping in mobile robotics.
He is post-doc researcher at the Geodetic Engineering Laboratory, École polytech-
nique fédérale de Lausanne (EPFL), Switzerland. His research interests lie in meth-
ods for high-accuracy navigation and mapping using small unmanned aerial vehi-
cles, but also span autonomous navigation and control, photogrammetry, computer
vision and signal processing.