A Novel 3D Reconstruction Sensor Using a Diving Lamp and a Camera for Underwater Cave Exploration

1. Introduction

Access to water has become one of the most pressing challenges facing humanity today. Yet water is abundant on a global scale: the oceans are an inexhaustible reserve, provided that salt and water are separated. Unfortunately, desalination methods are still too energy-intensive to be used profitably on a large scale. Pending a possible scientific breakthrough, we are therefore relying essentially on freshwater resources, which account for 2.5% of the world’s water [1] (Figure 1).

Only 1.2% of this freshwater is accessible at the surface, although it is now very polluted. Of the remainder, 68.7% is trapped in glaciers and ice caps and 30.1% is found underground. This groundwater, which represents 0.76% of the world’s resources, is extremely valuable. Groundwater is the result of a long and natural process of filtration through the various layers of soil, which means that its quality is generally considered to be very good. This groundwater is stored in aquifers, most of which are karstic.

Aquifer karstic structures, due to their complex and porous nature, present significant challenges in accurately mapping their intricate features. Traditional methods often rely on invasive techniques or sophisticated equipment, limiting accessibility and feasibility. In this paper, we propose a novel, non-invasive approach utilizing structured light projection in conjunction with a camera and a diving light to effectively map aquifer karstic formations. First part presents the main methods in underwater mapping. Second part details materials and methods of our original approach based on the principles of structured light, leveraging the precise projection of light contours onto the karstic surfaces. By capturing the resultant light patterns with a camera, we reconstruct detailed three-dimensional representations of the subsurface structures. Third part presents results to validate our approach by simulation of an aquifer galley first and finally through extensive outdoor experiments conducted in a particular equivalent aquifer karstic environment : The results demonstrate the efficacy of our method in accurately delineating intricate karstic features with remarkable detail and resolution. Furthermore, the non-destructive nature of this technique minimizes disturbance to delicate aquatic ecosystems while providing valuable insights into the subterranean landscape. This innovative methodology not only offers a cost-effective and non-invasive means of mapping aquifer karstic structures but also opens avenues for comprehensive environmental monitoring and resource management. Its potential applications span hydrogeological studies, environmental conservation efforts, and sustainable water resource management in karstic terrains worldwide.

2. Main Methods

Underwater mapping is a valuable tool in a wide range of applications, including marine biology, geology, archaeology, and the offshore industry. Various methods and sensors are employed to obtain 3D reconstructions of the environment, with data typically collected by divers, ROVs, or AUVs. In most cases, this involves performing surveys with downward-facing sensors and reconstructing the seafloor surface. The studies conducted by [2,3] provide comprehensive reviews of the various underwater 3D reconstruction methods. This synthesis of methods is illustrated in Figure 2.

Sensors used for 3D reconstruction can be classified into two main categories:

• active sensors: these sensors physically interact with the environment by emitting a signal or wave, such as sound waves in the case of sonars or light waves in the case of lidars.
• passive sensors: these sensors do not modify the environment, and the most commonly used ones for 3D reconstruction are camera-based.

Time-of-flight methods are the most direct methods for 3D reconstruction. They work on the same principle as a rangefinder: emitting a signal or wave and retrieving its echo. Knowing the speed of the signal thanks to the nature of the wave and the propagation medium, it is possible to proportionally deduce the relative distance between the sensor and the observed element from the measurement of the time elapsed between the emission of the wave and the return of the echo.

Among time-of-flight methods, those based on acoustic waves are by far the most used because their propagation is favored by the medium. They offer an interesting coverage capacity in front of the size of the areas to be studied in underwater mapping.

In addition and to allow the localization of 3D data during the measurement, most underwater navigation algorithms are based on acoustic positioning systems such as Doppler Velocity Log (DVL) and Ultra-short baseline (USBL) [4,5]. Multibeam sonars (MBS) are on the other hand adapted to bathymetric mapping (measurement of the seabed relief) [6] but much less to structured (non-flat) environments such as karst aquifers.

[7,8] presents a 3D model of the Zacaton cenote (a water-filled cave at least 300 meters deep) using sonar data collected by the DEPTHX (Deep Phreatic THermal eXplorer) vehicle. It is a NASA autonomous underwater vehicle measuring just over 2m and weighing over a ton, prefiguring among the robots that could explore underwater environments on other planets. Its sensors consist of a DVL, a sonar, an IMU and its movements are ensured by six motors allowing it to move in all directions.

Another method using acoustics is described by [9], which shows the results of an AUV performing limited penetration, inside an underwater cave using a mechanically scanned imaging sonar (MSIS for mechanically scanned imaging sonar). The results, obtained only for localization, show that it is possible to use this type of technology in our karst environment.

Efficient for obtaining bathymetry and macro-characteristics of habitats, these sensors nevertheless provide little information on fine characteristics (micro-bathymetry), do not allow to recover the texture of the environment and are very expensive.

Still in the class of time-of-flight methods, those based on electromagnetic waves, such as LiDAR, offer access to a very dense reconstruction of the seabed. However, on the other hand, they are reserved for very local observations, because the propagation of these waves are limited by the medium as we explained in the constraints. They are therefore not adapted to our situation.

The other class of methods is the triangulation method that we favor. These methods require at least two devices (or 2 different views) which will provide distance information from the same point, and thus form a triangle between the 2 devices and the point to be measured.

There are several ways to perform this triangulation, including:

• Multi-view triangulation involves taking multiple views of the same scene. By knowing the disparities of each viewpoint, it is possible to project the points of the scene onto different "image" planes, and thus triangulate the points of the scene to obtain a depth map. This technique is based on epipolar geometry explained in [10] with two possible approaches:

  -

  the use of a binocular stereo system where the relative positions of the two cameras are known. It is possible to make the system active, by using a structured light projector to add discriminant points or features related to the observed environment and thus facilitate the matching of points between the two images for triangulation.

  -

  the method known as Structure from Motion (SfM) involves taking a sequence of sequential images of an object or scene, typically from a single camera. To obtain a metric reconstruction, it is necessary to know the camera trajectory to resolve the scale ambiguity.

• structured light triangulation relies on an emitter-receiver system. The emitter is a light source (laser or not) that projects a pattern (distinctive patterns) onto a scene observed by the receiver, a camera. Using the positions of the camera and the projector, it is possible to triangulate the discriminant points brought by the light in the scene captured by the camera [11].

Regarding structured light triangulation methods, they are very accurate, but are generally only used to reconstruct small areas or objects. This is because light is quickly absorbed in water; therefore, we must be close to the target to best detect the projected light. Additionally, structured light sensors are not yet well suited for mobility and are therefore mostly used statically for the moment.

Cameras, on the other hand, have the advantage of being inexpensive, providing dense data as well as texture. However, most camera-based methods operate in open areas with natural lighting (or with artificial lighting that completely illuminates the workspace) and generally try to reconstruct the ground or small objects.

For example, [12,13] propose a 3D reconstruction method based on a stereo pair for archaeological objects in shallow water and in very poor conditions. Special filters are implemented to solve the noise problems caused by these poor conditions. This static 3D measurement is fused with a 3D map of the excavation area obtained using a scanning sonar. [14] performs a dense reconstruction of submerged structures at very high depths (up to 6000 meters) with a stereo system mounted on a controlled manipulator arm. The imposition of a known trajectory allows to know the different positions of the cameras to calculate a precise 3D model of the scene.

[15] reconstructs coral reefs with a wide-baseline stereo system using high-resolution video and a dense reconstruction method.

Our context is different from these works since we are in total darkness in a structured and extended environment whose content appears as the diver or robot moves.

[16] presents a description of the problems related to metrology in the underwater environment. The importance of calibrating the stereo pair underwater is emphasized, taking into account the refraction problems related to the medium. The detection of interest points is done on sufficiently bright images allowing a robust matching then used for a 3D reconstruction. The comparison of the object measurements obtained is made with a reference metrology and gives good results but in statics.

The work of [17] is directly related to our theme. The objective is to reconstruct a karstic underwater tunnel of Woodville Karst Plain in Florida using a new approach to stereovision generating a 3D point cloud of the cave. In their method, after calibration and rectification of the highly distorted images, they extract the contours of the artificial light projected into the cave which creates a cone of light illuminating part of the walls. These contours are then used as inputs in their stereovision algorithm. The displacement is estimated by visual odometry (ORB-SLAM) and allows a 3D reconstruction of a portion of the cave of about one hundred meters but without specifying the comparison of the result with a known ground truth.

For an even more optimal mapping of this underwater gallery, [18] have developed a SLAM method where they fuse this stereo method with a sonar and an inertial measurement unit, thus improving the accuracy of their measurements.

[19] uses a SFM method coupled with multi-view stereo photography to perform measurements of silting areas at the outlet of water transport pipelines to the seabed. The accuracy of their method allows a robust estimation of these elements useful for monitoring these areas but also remains static measurement.

[20] provides an overview of calibration and underwater image processing methods in structured light. [21] proposes a method of 3D reconstruction in underwater environments combining a SLAM approach and dense reconstruction using a stereo system embedded on the Aqua2 robot whose IMU allows localization. The results are mainly on coral reefs and partially on the reconstruction of underwater cavities. The question of the accuracy of the measurements with respect to a ground truth remains. [22] proposes an application of 3D reconstruction for underwater archeology. The system consists of two monochrome cameras and a sinusoidal light pattern projector. A combination of 3D reconstruction in structured light and SLAM-based visual odometry is proposed. All mounted on an ROV whose IMU data complements that of visual odometry. The first results are obtained on objects in a controlled environment. A test is then carried out over a dozen meters to follow a pseudo pipeline. But the size of the system is too large for an application in underground galleries.

[23] uses Ariane’s thread to enable the movement and guidance of a BlueRov robot. This learning method gives good results for localisation and assisted piloting. The next steps envisaged are autonomous navigation and then mapping of the cavities studied (but already explored by divers).

[24] proposes a real-time monocular 3D reconstruction method for underwater environments. It relies on optical flow tracking of characteristic seabed points and uses Delaunay triangulation to complete the seabed estimation. The method assumes a nearly planar seabed but remains unclear about the achieved metric accuracy.

3. Materials and Methods

The complete development of our method is detailed in Figure 3 and is based on three main blocks :

• calibration of the camera is used to estimate its intrinsic parameters, including radial distortion,
• calibration of the projector consists of estimating the geometric parameters of the cone (vertex ${\mathit{O}}_{\mathcal{P}}$, direction $\mathit{d}$, half-angle of opening $\alpha$),
• local 3D reconstruction based on the camera / projector triangulation leading to a 3D point cloud expressed in $\mathcal{C}$, the camera frame.

This triangulation is illustrated in Figure 4. It shows our system, with one camera and light projector. To simplify the figure, the light is projected onto a plane, but you can imagine that the procedure has to be applied to non-planar surfaces such as the walls of a gallery. Once the light contour has been detected in the image, the associated camera ray is calculated for each point on the contour. The Figure 4 shows two points ${\mathit{X}}^{\prime }$ and ${\mathit{Y}}^{\prime }$ belonging to the contour detected in the image and the rays associated with them. In the configuration shown in the Figure 4, the ray associated with ${\mathit{X}}^{\prime }$ (respectively ${\mathit{Y}}^{\prime }$) intersects the cone at two points ${\mathit{X}}_{\mathbf{1}}$ and ${\mathit{X}}_{\mathbf{2}}$ (respectively ${\mathit{Y}}_{\mathbf{1}}$ and ${\mathit{Y}}_{\mathbf{2}}$). To calculate these intersections we need to know the angle of aperture $\alpha$ and the transformation between the references of the cone and the camera called $T_{\{\mathcal{P}\to \mathcal{C}\}}$.

Still on the Figure 4, the intersections which correspond to the 3D points of the contour (${\mathit{X}}_{\mathbf{1}}$ and ${\mathit{Y}}_{\mathbf{2}}$ ) are shown in green and those which should not be taken into account (${\mathit{Y}}_{\mathbf{1}}$) are shown in red. The point ${\mathit{X}}_{\mathbf{2}}$ is not shown, but we can imagine that it would be red and beyond the area of the plane in the Figure 4. For a given ray, it is therefore necessary to identify whether it is the first intersection or the second intersection that belongs to the observed light contour.

This method of selecting intersections requires a study of the geometric relationships between the cone and the camera.

3.1. Cone Parameters

3.1.1. Cone of Revolution

The cone of revolution, referred to as a cone in the remainder of this document, is associated with its surface of revolution and not with its solid, as may be the case in certain applications. The terms we will use to define the relationship between a point and the cone are as follows:

• a point belonging to the cone is a point which belongs to the surface of the cone,
• a point inside the cone is a point which belongs to the solid bounded by the surface of the cone,
• a point outside the cone is a point which does not belong to the solid bounded by the surface of the cone.

This surface is generated by the revolution of a line, called the generatrix, which we will note g and which passes through a vertex, in this case ${\mathit{O}}_{\mathcal{P}}$.

The revolution takes place around a fixed axis which also passes through the vertex and which turns out to be an axis of symmetry of the cone, the direction of which is defined by the unit vector $\mathit{d}$.

If a plane not passing through ${\mathit{O}}_{\mathcal{P}}$ is orthogonal to the axis of symmetry, then its intersection with the cone is a circle.

Thus, the angle formed by g and the axis of symmetry is constant and corresponds to the half-angle of the cone opening $\alpha \in ]0;\frac{\pi }{2}[$.

The cone is represented in Figure 5 with the orthonormal base $({\mathit{O}}_{\mathcal{P}},{\mathit{x}}_{\mathcal{P}},{\mathit{y}}_{\mathcal{P}},{\mathit{z}}_{\mathcal{P}})$, oriented such that ${\mathit{z}}_{\mathcal{P}}=\mathit{d}$, which defines the $\mathcal{P}$ reference frame associated with the cone.

The generatrix g has the unit direction vector $\mathit{w}$ which is set by a new angle $\theta \in [0;2\pi [$ (Figure 5).

In the cone frame $\mathcal{P}$, we can simply express the vector ${\mathit{w}}_{\left\{\mathcal{P}\right\}}$ which then depends only on the angles $\theta$ and $\alpha$ :

$${\mathit{w}}_{\left\{\mathcal{P}\right\}}=\left[\begin{array}{c}sin\left(\alpha \right)cos\left(\theta \right)\\ sin\left(\alpha \right)sin\left(\theta \right)\\ cos\left(\alpha \right)\end{array}\right]$$

(1)

This vector can be used to express all the points on the surface parametrically.

Let ${\mathit{X}}_{\left\{\mathcal{P}\right\}}={\left[\begin{array}{ccc}{x}_{{\mathit{X}}_{\left\{\mathcal{P}\right\}}}& y_{{\mathit{X}}_{\left\{\mathcal{P}\right\}}}& z_{{\mathit{X}}_{\left\{\mathcal{P}\right\}}}\end{array}\right]}^T$ be a point $\mathit{X}$ expressed in the reference frame $\mathcal{P}$. Then, if $\mathit{X}$ belongs to the cone we have :

$${\mathit{X}}_{\left\{\mathcal{P}\right\}}=l{\mathit{w}}_{\left\{\mathcal{P}\right\}}\iff \left\{\begin{array}{c}\begin{array}{cc}\hfill & x_{{\mathit{X}}_{\left\{\mathcal{P}\right\}}}=lsin\left(\alpha \right)cos\left(\theta \right)\hfill \\ \hfill & y_{{\mathit{X}}_{\left\{\mathcal{P}\right\}}}=lsin\left(\alpha \right)sin\left(\theta \right)\hfill \\ \hfill & z_{{\mathit{X}}_{\left\{\mathcal{P}\right\}}}=lcos\left(\alpha \right)\hfill \end{array}\hfill \end{array}\right.$$

(2)

with l a real parameter whose absolute value defines the distance between ${\mathit{X}}_{\left\{\mathcal{P}\right\}}$ and the cone vertex. However, with $l\in \mathbb{R}$, the parametric equation defines a "double cone". However, since the cone models our light projector, we are only interested in the upper part of the "double cone", i.e. when the parameter $l\in {\mathbb{R}}^{+}$.

3.1.2. Orthogonal Distance

This parametrisation of the cone will help us to write the equation for the orthogonal distance between a point and the cone, an equation that will be useful for cone estimation.

Let $\mathit{M}$ be a point outside the cone and $\mathit{H}$ its orthogonal projection on the cone (see Figure 5).

The angle ${\theta }_H$ is the angle that gives the generatrix of the cone closest to $\mathit{M}$. This generatrix, called $g_H$, passes through $\mathit{H}$ and has the unit direction vector ${\mathit{w}}_H$ whose expression in the $\mathcal{P}$ reference frame is :

$${{\mathit{w}}_H}_{\left\{\mathcal{P}\right\}}=\left[\begin{array}{c}sin\left(\alpha \right)cos\left({\theta }_H\right)\\ sin\left(\alpha \right)sin\left({\theta }_H\right)\\ cos\left(\alpha \right)\end{array}\right]$$

If we also express $\mathit{M}$ in $\mathcal{P}$ with ${\mathit{M}}_{\left\{\mathcal{P}\right\}}={\left[\begin{array}{ccc}{x}_{{\mathit{M}}_{\left\{\mathcal{P}\right\}}}& y_{{\mathit{M}}_{\left\{\mathcal{P}\right\}}}& z_{{\mathit{M}}_{\left\{\mathcal{P}\right\}}}\end{array}\right]}^T$, we obtain the following expressions:

$$\left\{\begin{array}{c}\begin{array}{cc}\hfill {\theta }_H& =arctan\left(\frac{y_{{\mathit{M}}_{\left\{\mathcal{P}\right\}}}}{x_{{\mathit{M}}_{\left\{\mathcal{P}\right\}}}}\right)\hfill \\ \hfill {\mathit{H}}_{\left\{\mathcal{P}\right\}}& =\left({{\mathit{M}}_{\left\{\mathcal{P}\right\}}}^T·{{\mathit{w}}_H}_{\left\{\mathcal{P}\right\}}\right){{\mathit{w}}_H}_{\left\{\mathcal{P}\right\}}\hfill \end{array}\hfill \end{array}\right.$$

(3)

If ${{\mathit{M}}_{\left\{\mathcal{P}\right\}}}^T·{{\mathit{w}}_H}_{\left\{\mathcal{P}\right\}}<0$ then the projected point is in the lower part of the cone. As we are only working with the upper part of the cone, when this condition is true, the nearest point is the vertex of the ${\mathit{O}}_{\mathcal{P}}$ cone (see example Figure 6). The orthogonal distance h is therefore :

$$h=\left\{\begin{array}{c}\begin{array}{cccc}\hfill & \sqrt{{({\mathit{M}}_{\left\{\mathcal{P}\right\}}-{\mathit{H}}_{\left\{\mathcal{P}\right\}})}^T({\mathit{M}}_{\left\{\mathcal{P}\right\}}-{\mathit{H}}_{\left\{\mathcal{P}\right\}})}\hfill & \hfill & \mathrm{si}{{\mathit{M}}_{\left\{\mathcal{P}\right\}}}^T{{\mathit{w}}_H}_{\left\{\mathcal{P}\right\}}⩾0\hfill \\ \hfill & \sqrt{{{\mathit{M}}_{\left\{\mathcal{P}\right\}}}^T{\mathit{M}}_{\left\{\mathcal{P}\right\}}}\hfill & \hfill & \mathrm{si}{{\mathit{M}}_{\left\{\mathcal{P}\right\}}}^T{{\mathit{w}}_H}_{\left\{\mathcal{P}\right\}}<0\hfill \end{array}\hfill \end{array}\right.$$

(4)

3.1.3. Quadratic form of the Cone

To obtain this relationship, we must first note that if the angle between the vectors $(\mathit{X}-{\mathit{O}}_{\mathcal{P}})$ and $\mathit{d}$ or $(\mathit{X}-{\mathit{O}}_{\mathcal{P}})$ and $-\mathit{d}$ is equal to the angle $\alpha$, then the point $\mathit{X}$ belongs to the cone. This is equivalent to the following scalar product:

$${(\mathit{X}-{\mathit{O}}_{\mathcal{P}})}^T\mathit{d}=\pm ∥\mathit{X}-{\mathit{O}}_{\mathcal{P}}∥cos\left(\alpha \right)$$

When squared, the expression becomes :

$$\begin{array}{cc}\hfill & {\left({(\mathit{X}-{\mathit{O}}_{\mathcal{P}})}^T\mathit{d}\right)}^2={\left(∥\mathit{X}-{\mathit{O}}_{\mathcal{P}}∥cos\left(\alpha \right)\right)}^2\hfill \\ \hfill \iff & \left({(\mathit{X}-{\mathit{O}}_{\mathcal{P}})}^T\mathit{d}\right)\left({(\mathit{X}-{\mathit{O}}_{\mathcal{P}})}^T\mathit{d}\right)={(\mathit{X}-{\mathit{O}}_{\mathcal{P}})}^T(\mathit{X}-{\mathit{O}}_{\mathcal{P}}){cos}^2\left(\alpha \right)\hfill \\ \hfill \iff & {(\mathit{X}-{\mathit{O}}_{\mathcal{P}})}^T\mathit{d}{\mathit{d}}^T(\mathit{X}-{\mathit{O}}_{\mathcal{P}})-{(\mathit{X}-{\mathit{O}}_{\mathcal{P}})}^T(\mathit{X}-{\mathit{O}}_{\mathcal{P}}){cos}^2\left(\alpha \right)=0\hfill \\ \hfill \iff & {(\mathit{X}-{\mathit{O}}_{\mathcal{P}})}^T\underset{Q}{\underbrace{\left[\begin{array}{c}\mathit{d}{\mathit{d}}^T-{cos}^2\left(\alpha \right)I_3\end{array}\right]}}(\mathit{X}-{\mathit{O}}_{\mathcal{P}})=0\hfill \end{array}$$

We thus obtain a quadratic function of the cone in the general case, directly showing the parameters $\mathit{d}$, $\alpha$ and ${\mathit{O}}_{\mathcal{P}}$ which define it:

$$\mathcal{Q}\left(\mathit{X}\right)={(\mathit{X}-{\mathit{O}}_{\mathcal{P}})}^{T}Q(\mathit{X}-{\mathit{O}}_{\mathcal{P}})$$

(5)

This gives the relationships between the point $\mathit{X}$ and the cone:

3.2. Projector Calibration

Projector calibration is an essential step, since it involves estimating the parameters of the cone used in our method and in the equations we have seen so far. As a reminder, the parameters of the cone to be estimated are its vertex ${\mathit{O}}_{\mathcal{P}}$, its direction vector $\mathit{d}$ and its half-angle of aperture $\alpha$.

Since ${\mathit{O}}_{\mathcal{P}}$ and $\mathit{d}$ represent the relative pose of the cone with respect to the camera, the projector must be fixed with respect to the camera during calibration and during any experimentation. Furthermore, if the projector has a variable half-angle aperture, it is also necessary to ensure that this is fixed.

3.2.1. 3D Point Generation for Cone Estimation

The first step in the calibration method is to generate a cloud of 3D points belonging to the cone. To achieve it, the chosen method is to capture several light projection images on a flat surface.

This light projection is the result of the intersection between the projector cone and the flat surface, which is by definition a conic. In our case, it will be an ellipse because we want a closed shape. The projection of this ellipse into the image is also an ellipse.

If we can estimate the relative pose of this surface with respect to the camera and extract the light contour in the image, which is an ellipse, then we can obtain the real ellipse from the ellipse in the image.

To estimate the relative pose, we attached a chessboard pattern to the surface to estimate its relative pose. This relative pose is obtained by estimating the homography between the calibration pattern and the image plane [25]. As the camera is calibrated, rotation and translation can be estimated from the homography.

To extract the ellipse from the image, we again use the contour detection method presented in [26,27].

Then, from the contour points obtained, we estimate the best ellipse using the method of [28] which is based on least squares optimisation.

It is therefore possible to obtain n elliptical sections of the cone from n images taken at different distances from the surface. For each ellipse extracted in an image, we obtain the ellipse in the surface via the estimated homography between the image plane and this surface. The estimated relative pose allows us to obtain the elliptical section, i.e. the set of 3D points expressed in the camera frame of reference belonging to the ellipse of the surface.

Figure 8 shows an example where the camera has captured three images of an ellipse surrounding a chessboard pattern on the flat surface in three different poses. This is equivalent to obtaining three elliptical sections of the same cone.

Figure 7. Example of detecting the contour of a halo of light on a wall using a chessboard to estimate the plane pose. The contour is in green and the adjusted ellipse in blue.

Figure 7. Example of detecting the contour of a halo of light on a wall using a chessboard to estimate the plane pose. The contour is in green and the adjusted ellipse in blue.

Consequently, if we have obtained n elliptical sections where m 3D points are extracted per section, we generate a cloud of $n\times m$ 3D points belonging to the cone.

3.2.2. Cone Estimation

The aim is now to estimate the cone that best approximates the 3D point cloud. This is done by geometric fitting, i.e. by minimising the orthogonal distances between the 3D points and the cone.

Let $\mathit{p}$ be the vector containing the six parameters of the cone, which are :

• the three coordinates of its vertex ${\mathit{O}}_{\mathcal{P}}$,
• the two angles yaw and pitch in ZYX-Euler convention of its direction vector $\mathit{d}$, the angle roll not being necessary since a cone has an axis of symmetry,
• its opening half-angle $\alpha$.

To estimate the $\mathit{p}$ vector, we use the Levenberg-Marquardt (LM) iterative method [29]. For the first iteration, we need to choose an initial parameter vector ${\mathit{p}}_0$. This could be the null vector, for example. In practice, we will use measurements taken with a caliper for the various parameters of ${\mathit{p}}_0$.

At each iteration i, LM receives a solution vector ${\mathit{p}}_i$ and the function for calculating the vector of orthogonal distances to be minimised. Our function first calculates the rotation and translation of the cone relative to the camera from ${\mathit{O}}_{\mathcal{P}}$ and the yaw and pitch angles of the vector ${\mathit{p}}_i$. It then expresses the 3D points in the cone reference frame and calculates the orthogonal distance vector between the 3D points and the cone using the equation 4.

So from ${\mathit{p}}_i$ and this function, LM returns a new solution vector ${\mathit{p}}_{i+1}$.

When the difference between the sum of the squared orthogonal distances obtained with ${\mathit{p}}_j$ and the calculus obtained with ${\mathit{p}}_{j+1}$ is less than an arbitrary threshold, this means that LM has reached convergence at iteration j. The best solution in the least squares sense is therefore the vector ${\mathit{p}}_j$.

It is important to note that there are an infinite number of cones that share the same elliptical section. On the other hand, there is only one cone that passes through an elliptical section and through a point outside that section. Consequently, this minimisation can only work if we have at least two elliptical sections.

3.3. Intersection with a Ray

Now that we know how to determine the cone parameters using our calibration method, the next step is to be able to calculate the intersections between a camera ray and the cone. We are only interested in the intersections with the top of the cone.

A ray is a half-line that can be expressed parametrically using the following equation:

$$\mathit{X}\left(t\right)=\mathit{O}+t\mathit{u}$$

(7)

with $\mathit{O}$ the starting point of the ray, $\mathit{u}$ its unit direction vector and t a positive real number such that $\mathit{X}\left(0\right)=\mathit{O}$. Note that we will only deal with cases where the point $\mathit{O}$ is outside the cone. This is because the starting point of the camera rays is the optical centre in our situation. The point $\mathit{O}$ can nevertheless be inside the lower part of the cone.

If $\mathit{X}\left(t\right)$ is an intersection point, then it belongs to the cone. Therefore, to determine the intersections, we can simply calculate the values of t which cancel the general polynomial of the cone (see 5) after replacing $\mathit{X}$ by the expression for $\mathit{X}\left(t\right)$ :

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill & \mathcal{Q}\left(\mathit{X}\left(t\right)\right)={(\mathit{X}\left(t\right)-{\mathit{O}}_{\mathcal{P}})}^{T}Q(\mathit{X}\left(t\right)-{\mathit{O}}_{\mathcal{P}})=0\hfill \\ \hfill & \iff {(\mathit{O}+t\mathit{u}-{\mathit{O}}_{\mathcal{P}})}^{T}Q(\mathit{O}+t\mathit{u}-{\mathit{O}}_{\mathcal{P}})=0\hfill \\ \hfill & \iff t^2\left({\mathit{u}}^{T}Q\mathit{u}\right)+t\left(2{\mathit{u}}^{T}Q(\mathit{O}-{\mathit{O}}_{\mathcal{P}})\right)+{(\mathit{O}-{\mathit{O}}_{\mathcal{P}})}^{T}Q(\mathit{O}-{\mathit{O}}_{\mathcal{P}})=0\hfill \\ \hfill & \iff a{t}^2+bt+c=0\mathrm{with}\left\{\begin{array}{c}\begin{array}{cc}\hfill a& ={\mathit{u}}^{T}Q\mathit{u}\hfill \\ \hfill b& =2{\mathit{u}}^{T}Q(\mathit{O}-{\mathit{O}}_{\mathcal{P}})\hfill \\ \hfill c& ={(\mathit{O}-{\mathit{O}}_{\mathcal{P}})}^{T}Q(\mathit{O}-{\mathit{O}}_{\mathcal{P}})=\mathcal{Q}\left(\mathit{O}\right)\hfill \end{array}\hfill \end{array}\right.\hfill \end{array}\end{array}$$

3.4. Projection of the Generatrices of the Cone in the Image

To understand the relationship between a point on the cone and its projection into the image, we need to look at the projections of the generatrices of the cone. Since the cone represents the projector, we will only consider the upper half of the cone, which means that the generatrices can be limited to the half-lines that start at the vertex ${\mathit{O}}_{\mathcal{P}}$. Note that the study is limited to the case where the cone is oriented in the same direction as the camera, as required by our method. Geometrically, this condition is satisfied if the angle called $\mu$ between ${\mathit{z}}_{\mathcal{C}}$ and $\mathit{d}={\mathit{z}}_{\mathcal{P}}$ is in the interval ]$\frac{\pi }{2}-\alpha ;-\frac{\pi }{2}+\alpha$[ (see Figure 9). In practice, for the camera to be able to see the light projection, the angle $\mu$ will be close to 0.

3.4.1. Calculating the Projections of the Generatrices of the Cone

Let us call ${\Pi }_g$ the plane which passes through the optical centre ${\mathit{O}}_{\mathcal{C}}$ and through at least one generatrix of the cone. The intersection of this plane with the image plane is the line ${\ell }_g$.

Let ${\mathit{X}}^{\prime }$ be a point in the image such that ${\mathit{X}}^{\prime }=K{\mathcal{C}}_{\left\{\mathit{u}\right\}}$ with $\mathit{u}$ the direction vector of the camera ray associated with ${\mathit{X}}^{\prime }$. If ${\mathit{X}}^{\prime }$ is the projection into the image of a 3D point belonging to a generatrix resulting from the intersection between the cone and the plane ${\Pi }_g$, then ${\mathit{X}}^{\prime }$ belongs to the line ${\ell }_g$. However, if ${\mathit{X}}^{\prime }$ belongs to the line ${\ell }_g$, it is not necessarily the projection of a point on the cone. To obtain the equation of ${\ell }_g$, the simplest way is to select two points belonging to one of the generatrices (or the generatrix) through which the ${\Pi }_g$ passes, then project them into the image to find the parameters of the line. The same result would be obtained by projecting two points of the plane ${\Pi }_g$ which do not belong to the same camera radius. However, knowing the line ${\ell }_g$ is not enough to understand the projection of a generatrix.

Observation When ${\Pi }_g$ Is Not Tangent to the Cone

When the plane ${\Pi }_g$ is not tangent to the cone, it passes through two generatrices $g_1$ and $g_2$, with $g_1$ chosen as the generatrix closest to the optical centre ${\mathit{O}}_{\mathcal{C}}$ (see Figure 10).

The direction vectors of these two generatrices are ${\mathit{w}}_{\mathbf{1}}$ and ${\mathit{w}}_{\mathbf{2}}$ such that :

$${{\mathit{w}}_{\mathbf{1}}}_{\left\{\mathcal{P}\right\}}=\left[\begin{array}{c}sin\left(\alpha \right)cos\left({\theta }_1\right)\\ sin\left(\alpha \right)sin\left({\theta }_1\right)\\ cos\left(\alpha \right)\end{array}\right]\mathrm{et}{{\mathit{w}}_{\mathbf{2}}}_{\left\{\mathcal{P}\right\}}=\left[\begin{array}{c}sin\left(\alpha \right)cos\left({\theta }_2\right)\\ sin\left(\alpha \right)sin\left({\theta }_2\right)\\ cos\left(\alpha \right)\end{array}\right]$$

with ${\theta }_1$ (respectively ${\theta }_2$) the angle between the vector ${\mathit{w}}_{\mathbf{1}}$ (respectively ${\mathit{w}}_{\mathbf{2}}$) projected in the plane ${\mathit{x}}_{\mathcal{P}}{\mathit{y}}_{\mathcal{P}}$ and the vector ${\mathit{x}}_{\mathcal{P}}$ (see Figure 10). As the plane ${\Pi }_g$ is chosen here as not tangent to the cone, we have ${\theta }_1\ne {\theta }_2$. In this configuration, if ${\mathit{X}}^{\prime }$ is a point on the line ${\ell }_g$ then it can be both the projection of a point ${\mathit{X}}_{\mathbf{1}}$ belonging to $g_1$ and of a point ${\mathit{X}}_{\mathbf{2}}$ belonging to $g_2$. The generatrix $g_1$ therefore contains the first intersections of the camera ray with the cone, while $g_2$ contains the second intersections.

Calculating the Two Special Generatrices When ${\Pi }_g$ Is Tangent to the Cone

There are only two orientations for which the plane ${\Pi }_g$ is tangent to the cone. When this happens, then ${\theta }_1={\theta }_2$ and the generatrices $g_1$ and $g_2$ are superimposed. The two possible cases are illustrated in Figure 11. The two angles ${\theta }_A$ and ${\theta }_B$ define the two generatrices $g_A$ and $g_B$ through which pass the two planes tangent to the cone ${\Pi }_{g_A}$ and ${\Pi }_{g_B}$. It is important to note that the angles ${\theta }_A$ and ${\theta }_B$ depend solely on the parameters of the cone and its pose relative to the camera $T_{\{\mathcal{E}\to \mathcal{C}\}}$, which means that they are unique and fixed for a given configuration.

To determine these two angles, we will use a plane tangent to the cone called ${\Pi }_T$ which passes through a generatrix g whose direction vector is $\mathit{w}$ defined by an angle $\theta$. It is represented in Figure 10 passing through the generatrix $g_2$ purely for reasons of visibility. The orientation of the plane ${\Pi }_T$ is characterised by the unit normal vector ${\mathit{n}}_{\mathit{T}}$ such that :

$${{\mathit{n}}_{\mathit{T}}}_{\left\{\mathcal{P}\right\}}=\left[\begin{array}{c}cos\left(\alpha \right)cos\left(\theta \right)\\ cos\left(\alpha \right)sin\left(\theta \right)\\ -sin\left(\alpha \right)\end{array}\right]$$

(8)

We are therefore looking for the two solutions of $\theta$ so that the plane ${\Pi }_T$ is superimposed on the plane ${\Pi }_{g_A}$ or ${\Pi }_{g_B}$. This superposition only occurs when the plane ${\Pi }_T$ passes through the optical centre ${\mathit{O}}_{\mathcal{C}}$ and therefore when the vector ${\mathit{n}}_{\mathit{T}}$ is orthogonal to the vector (${\mathit{O}}_{\mathcal{C}}$ - ${\mathit{O}}_{\mathcal{P}}$). We are therefore looking for $\theta$ when :

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill & {{{\mathit{O}}_{\mathcal{C}}}_{\left\{\mathcal{P}\right\}}}^T{{\mathit{n}}_{\mathit{T}}}_{\left\{\mathcal{P}\right\}}=0\iff \left[\begin{array}{ccc}{x}_{{{\mathit{O}}_{\mathcal{C}}}_{\left\{\mathcal{P}\right\}}}& y_{{{\mathit{O}}_{\mathcal{C}}}_{\left\{\mathcal{P}\right\}}}& z_{{{\mathit{O}}_{\mathcal{C}}}_{\left\{\mathcal{P}\right\}}}\end{array}\right]\left[\begin{array}{c}cos\left(\alpha \right)cos\left(\theta \right)\\ cos\left(\alpha \right)sin\left(\theta \right)\\ -sin\left(\alpha \right)\end{array}\right]=0\hfill \\ \hfill \iff & x_{{{\mathit{O}}_{\mathcal{C}}}_{\left\{\mathcal{P}\right\}}}cos\left(\alpha \right)cos\left(\theta \right)+y_{{{\mathit{O}}_{\mathcal{C}}}_{\left\{\mathcal{P}\right\}}}cos\left(\alpha \right)sin\left(\theta \right)-z_{{{\mathit{O}}_{\mathcal{C}}}_{\left\{\mathcal{P}\right\}}}sin\left(\alpha \right)=0\hfill \end{array}\end{array}$$

(9)

with ${{\mathit{O}}_{\mathcal{C}}}_{\left\{\mathcal{P}\right\}}$ the optical centre expressed in the reference frame of the $\mathcal{P}$ projector. To solve the equation 9 we pose :

$$\begin{array}{ccccc}\hfill \mathit{s}=cos\left(\alpha \right)\left[\begin{array}{c}{x}_{{{\mathit{O}}_{\mathcal{C}}}_{\left\{\mathcal{P}\right\}}}\\ y_{{{\mathit{O}}_{\mathcal{C}}}_{\left\{\mathcal{P}\right\}}}\end{array}\right]& \hfill & \hfill \mathrm{and}& \hfill & \hfill \mathit{t}=\left[\begin{array}{c}cos\left(\theta \right)\\ sin\left(\theta \right)\end{array}\right]\end{array}$$

(10)

We can rewrite the equation 9 using the scalar product between the vectors $\mathit{s}$ and $\mathit{t}$ and the angles defined in Figure 12 :

$${\mathit{s}}^T\mathit{t}=\stackrel{∥\mathit{s}∥∥\mathit{t}∥=∥\mathit{s}∥}{\overbrace{cos\left(\alpha \right)\sqrt{x_{{{\mathit{O}}_{\mathcal{C}}}_{\left\{\mathcal{P}\right\}}}^2+y_{{{\mathit{O}}_{\mathcal{C}}}_{\left\{\mathcal{P}\right\}}}^2}}}cos\left(\gamma \right)=z_{{{\mathit{O}}_{\mathcal{C}}}_{\left\{\mathcal{P}\right\}}}sin\left(\alpha \right)$$

The angle $\gamma$ between the vectors $\mathit{s}$ and $\mathit{t}$ is therefore :

$$\gamma =\pm arccos\left(\frac{z_{{{\mathit{O}}_{\mathcal{C}}}_{\left\{\mathcal{P}\right\}}}sin\left(\alpha \right)}{cos\left(\alpha \right)\sqrt{x_{{{\mathit{O}}_{\mathcal{C}}}_{\left\{\mathcal{P}\right\}}}^2+y_{{{\mathit{O}}_{\mathcal{C}}}_{\left\{\mathcal{P}\right\}}}^2}}\right)=\pm arccos(\frac{z_{{{\mathit{O}}_{\mathcal{C}}}_{\left\{\mathcal{P}\right\}}}}{\sqrt{x_{{{\mathit{O}}_{\mathcal{C}}}_{\left\{\mathcal{P}\right\}}}^2+y_{{{\mathit{O}}_{\mathcal{C}}}_{\left\{\mathcal{P}\right\}}}^2}}tan\left(\alpha \right))$$

As for the angle $\beta$ which defines the orientation of the vector $\mathit{s}$, its expression is :

$$\beta =\mathrm{atan}2(y_{{{\mathit{O}}_{\mathcal{C}}}_{\left\{\mathcal{P}\right\}}},x_{{{\mathit{O}}_{\mathcal{C}}}_{\left\{\mathcal{P}\right\}}})$$

The solution to equation 9 is therefore :

$${\theta }_{A,B}=\beta \pm \gamma =\mathrm{atan}2(y_{{{\mathit{O}}_{\mathcal{C}}}_{\left\{\mathcal{P}\right\}}},x_{{{\mathit{O}}_{\mathcal{C}}}_{\left\{\mathcal{P}\right\}}})\pm arccos(\frac{z_{{{\mathit{O}}_{\mathcal{C}}}_{\left\{\mathcal{P}\right\}}}}{\sqrt{x_{{{\mathit{O}}_{\mathcal{C}}}_{\left\{\mathcal{P}\right\}}}^2+y_{{{\mathit{O}}_{\mathcal{C}}}_{\left\{\mathcal{P}\right\}}}^2}}tan\left(\alpha \right))$$

(11)

Note that the two angles ${\theta }_A$ and ${\theta }_B$ only depend on the half-angle of the opening $\alpha$ of the cone and the translation between the cone and the camera (translation expressed in the equation by ${{\mathit{O}}_{\mathcal{C}}}_{\left\{\mathcal{P}\right\}}$). From the expression of these two angles, we can obtain the expression of the two unit direction vectors ${\mathit{w}}_{\mathit{A}}$ and ${\mathit{w}}_{\mathit{B}}$ of the two special generatrices $g_A$ and $g_B$. The projections of these two generatrices belong to the lines ${\ell }_{g_A}$ and ${\ell }_{g_B}$ which are the intersections of the image plane with the planes ${\Pi }_{g_A}$ and ${\Pi }_{g_B}$.

Splitting the Cone into Two Areas Using the Two Special Generatrices

Using the generatrices $g_A$ and $g_B$ we can divide the cone into two surfaces as shown in Figure 13. The cyan surface named ${\mathcal{S}}_1$ (resp. magenta surface named ${\mathcal{S}}_2$) contains the set of generatrices containing all the first (resp. second) possible intersections between a camera beam and the cone, generatrices previously named $g_1$ (resp. $g_2$).

In the image plane, this shows that we can delimit the areas where the camera rays can intersect the cone. The cyan area, surrounded by the ellipse ${\mathcal{E}}_{\infty }^{\prime }$, contains the projections of part of the set of 3D points belonging to the surface ${\mathcal{S}}_1$. The magenta area contains the projections of the other part of the set of 3D points belonging to the surface ${\mathcal{S}}_1$ and contains the projections of the set of 3D points belonging to the surface ${\mathcal{S}}_2$. These areas are bounded by :

• the lines ${\ell }_{g_A}$ and ${\ell }_{g_B}$.
• the ellipse ${\mathcal{E}}_{\infty }^{\prime }$ which is the set of vanishing points of the generatrices. These points define a first bound on all the projections of the generatrices.
• the projection of the vertex ${\mathit{O}}_{\mathcal{P}}$ into the image plane, named ${\mathit{O}}_{\mathcal{P}}^{\prime }$. It defines a second bound on all the projections of the generatrices.

In the configuration illustrated in Figure 13, each generatrix projection is therefore a segment whose limits are the vanishing point of the generatrix belonging to the ellipse ${\mathcal{E}}_{\infty }^{\prime }$ and the point ${\mathit{O}}_{\mathcal{P}}^{\prime }$.

However, there are other configurations where the projections of the generatrices onto the image plane are not segments but lines. Indeed, if $z_{{{\mathit{O}}_{\mathcal{P}}}_{\left\{\mathcal{C}\right\}}}=0$, then the line passing through ${\mathit{O}}_{\mathcal{P}}$ and ${\mathit{O}}_{\mathcal{C}}$ is parallel to the plane, which implies that the projection of the vertex ${\mathit{O}}_{\mathcal{P}}$ into the image plane is a point at infinity (or ideal point) and that the lines ${\ell }_{g_A}$ and ${\ell }_{g_B}$ are parallel.

The case where $z_{{{\mathit{O}}_{\mathcal{P}}}_{\left\{\mathcal{C}\right\}}}<0$ is illustrated in Figure 14 which shows the lower part of the cone behind the camera. This part is shown in brown and its projection in the image plane is also in brown. This projection is obviously fictitious, as this part cannot be seen by the camera even if it had an infinite field of view. The point ${\mathit{O}}_{\mathcal{P}}^{\prime }$, also fictitious, is still the intersection of the lines ${\ell }_{g_A}$ and ${\ell }_{g_B}$, but this time it is on the other side in comparaison of the Figure 13.

As for the magenta and cyan areas in the image plane, they are still delimited by the lines ${\ell }_{g_A}$ and ${\ell }_{g_B}$ and the ellipse ${\mathcal{E}}_{\infty }^{\prime }$ but extend towards infinity on the side where the lines ${\ell }_{g_A}$ and ${\ell }_{g_B}$ diverge.

So if $z_{{{\mathit{O}}_{\mathcal{P}}}_{\left\{\mathcal{C}\right\}}}\le 0$ then each generatrix projection is a half-line whose only boundary is its vanishing point belonging to the ellipse ${\mathcal{E}}_{\infty }^{\prime }$. However, as we are working with real cameras and not with the image plane, which is infinite, each generatrix will have a segment as its projection into the image.

This division of the cone into two surfaces leads to several results:

• Any ray associated with a point in the cyan area has as its unique intersection with the cone a 3D point belonging to the surface ${\mathcal{S}}_1$.
• Any ray associated with a point in the magenta area has two intersections with the cone, the first of which belongs to the surface ${\mathcal{S}}_1$ and the second to the surface ${\mathcal{S}}_2$. These two intersections are superimposed if the point in the magenta area belongs to the line ${\ell }_{g_A}$ or to the line ${\ell }_{g_B}$.
• Any ray associated with a point outside the magenta area and the cyan area does not have any intersection with the cone.

3.5. Intersection Selection and Triangulation

We now have everything we need to develop a method for determining which intersection to choose when the ray associated with a pixel in the contour has two intersections with the cone. To better understand the various stages of the procedure, let’s take as an example the 3D contour of light whose associated 3D curve is circular and results from the intersection between a plane and the cone.

The example of the circular contour is illustrated in Figure 13 which contains the following elements:

• ${\mathcal{C}}^{\prime }$ is the closed 2D curve corresponding to the contour in the image with ${\mathit{X}}^{\prime }$ a 2D point on the curve.
• $\mathcal{C}$ is our circular 3D curve corresponding to the contour of the light in the scene, with $\mathit{X}$ a 3D point on the curve which has ${\mathit{X}}^{\prime }$ as its projection. The point $\mathit{X}$ therefore corresponds to the correct intersection between the cone and the camera ray associated with ${\mathit{X}}^{\prime }$.
• The areas containing the first and second intersections with the cone are delimited by the two generatrices $g_A$ and $g_B$ which divide the cone into two surfaces ${\mathcal{S}}_1$ and ${\mathcal{S}}_2$. They also divide $\mathcal{C}$ into two curves, one in cyan called ${\mathcal{C}}_1$ and the other in magenta called ${\mathcal{C}}_2$.
• $\mathit{A}$ and $\mathit{B}$ are respectively the two intersections of $g_A$ and $g_B$ with the curve $\mathcal{C}$ and are thus the only two 3D points common to the curves ${\mathcal{C}}_1$ and ${\mathcal{C}}_2$.
• The cyan curve named ${{\mathcal{C}}^{\prime }}_1$ and the magenta curve named ${{\mathcal{C}}^{\prime }}_2$ are respectively the projections of the curves ${\mathcal{C}}_1$ and ${\mathcal{C}}_2$.
• ${\mathit{A}}^{\prime }$ and ${\mathit{B}}^{\prime }$ are the projections of $\mathit{A}$ and $\mathit{B}$ and are therefore the only two points common to the curves ${{\mathcal{C}}^{\prime }}_1$ and ${{\mathcal{C}}^{\prime }}_2$. They belong to both the curve ${\mathcal{C}}^{\prime }$ and the lines ${\ell }_{g_A}$ and ${\ell }_{g_B}$. So, in theory, they correspond to the two unique points of tangency of the lines ${\ell }_{g_A}$ and ${\ell }_{g_B}$ with ${\mathcal{C}}^{\prime }$.

So, if we get the points of tangency ${\mathit{A}}^{\prime }$ and ${\mathit{B}}^{\prime }$, we will be able to divide the 2D curve ${\mathcal{C}}^{\prime }$ in two via the line $\left({\mathit{A}}^{\prime }{\mathit{B}}^{\prime }\right)$. However, in practice :

• the curve ${\mathcal{C}}^{\prime }$ is discrete,
• the light contours are extracted in a perfectible way,
• the camera calibration and the cone estimation have uncertainties, so the lines ${\ell }_{g_A}$ and ${\ell }_{g_B}$ also have uncertainties.

All this means that the lines ${\ell }_{g_A}$ and ${\ell }_{g_B}$ are not exactly tangent to the curve ${\mathcal{C}}^{\prime }$. The method for determining the points ${\mathit{A}}^{\prime }$ and ${\mathit{B}}^{\prime }$ must therefore be adapted to this constraint.

Once these two points have been obtained and the curve ${\mathcal{C}}^{\prime }$ has been divided in two, it remains to determine which of the two curves is ${{\mathcal{C}}^{\prime }}_1$ or ${{\mathcal{C}}^{\prime }}_2$. To do this, we need to look at the relative position of the vertex ${\mathit{O}}_{\mathcal{P}}$ in relation to the camera. Figure 15 shows that the vertex ${\mathit{O}}_{\mathcal{P}}$ and the curve ${{\mathcal{C}}^{\prime }}_2$ lie in the same half-space bounded by the plane passing through ${\mathit{O}}_{\mathcal{C}}$ and the line $\left({\mathit{A}}^{\prime }{\mathit{B}}^{\prime }\right)$. This property is always true whatever the pose of the cone.

In the image plane, this property can be expressed by defining a point ${\mathit{p}}_{\mathit{ref}}$ such that ${\mathit{p}}_{\mathit{ref}}$ and the curve ${{\mathcal{C}}^{\prime }}_2$ lie in the same half-plane bounded by the line $\left({\mathit{A}}^{\prime }{\mathit{B}}^{\prime }\right)$.

There are an infinite number of solutions for ${\mathit{p}}_{\mathit{ref}}$ in the image plane which satisfy this condition. We have therefore forced ${\mathit{p}}_{\mathit{ref}}$ to belong to the circumscribed circle of the rectangle delimiting the image, as can be seen in Figure 15. The segment between ${\mathit{p}}_{\mathit{ref}}$ and the centre of the image forms an angle with the horizontal axis of the image called $\nu$. This gives :

$${\mathit{p}}_{\mathit{ref}}=\sqrt{w^2+h^2}\left[\begin{array}{c}cos\left(\nu \right)\\ sin\left(\nu \right)\end{array}\right]$$

(12)

with w and h respectively the width and height of the image in pixels. For the condition to be true, one possible solution for $\nu$ is to indicate the direction of ${\mathit{O}}_{\mathcal{P}}$ relative to the optical centre ${\mathit{O}}_{\mathcal{C}}$. To obtain this angle, we orthogonally projected ${\mathit{O}}_{\mathcal{P}}$ onto the plane normal to the camera axis passing through ${\mathit{O}}_{\mathcal{C}}$, which gives the point ${\mathit{H}}_{\mathcal{P}}$ illustrated in Figure 15. This gives $\nu$ which is the angle formed by the vectors $({\mathit{H}}_{\mathcal{P}}-{\mathit{O}}_{\mathcal{C}})$ and ${\mathit{x}}_{\mathcal{C}}$ :

$$({\mathit{H}}_{\mathcal{P}}-{\mathit{O}}_{\mathcal{C}})=\left[\begin{array}{c}{x}_{{{\mathit{O}}_{\mathcal{P}}}_{\left\{\mathcal{C}\right\}}}\\ y_{{{\mathit{O}}_{\mathcal{P}}}_{\left\{\mathcal{C}\right\}}}\\ 0\end{array}\right]\mathrm{donc}\nu =arctan\left(\frac{y_{{{\mathit{O}}_{\mathcal{P}}}_{\left\{\mathcal{C}\right\}}}}{x_{{{\mathit{O}}_{\mathcal{P}}}_{\left\{\mathcal{C}\right\}}}}\right)$$

(13)

Calculating this angle $\nu$ and the point ${\mathit{p}}_{\mathit{ref}}$ allows us to deduce which of the two curves is ${{\mathcal{C}}^{\prime }}_2$.

Once the curves ${{\mathcal{C}}^{\prime }}_1$ and ${{\mathcal{C}}^{\prime }}_2$ have been obtained, it is then possible to obtain the 3D contour. Indeed, if a point ${\mathit{X}}^{\prime }$ belongs to ${{\mathcal{C}}^{\prime }}_1$, then the first intersection of the camera ray associated with it must be chosen, and if it belongs to ${{\mathcal{C}}^{\prime }}_2$, the second must be chosen. This is because any point ${\mathit{X}}^{\prime }$ of ${{\mathcal{C}}^{\prime }}_1$ is the image of a 3D point $\mathit{X}$ of ${\mathcal{C}}_1$, and any point ${\mathit{X}}^{\prime }$ of ${{\mathcal{C}}^{\prime }}_2$ is the image of a 3D point $\mathit{X}$ of ${\mathcal{C}}_2$.

3.6. Test of the Method in Simulation

Contour Simulation

To illustrate the different stages of the method, we will use a gallery model. Figure 16 shows the camera represented by a pyramid and the projector represented by a cone positioned inside our model. In our example, the axis of the cone is parallel to the axis of the camera. The vertex ${\mathit{O}}_{\mathcal{P}}$ is positioned above and to the right of the camera ; it belongs to the plane which passes through the optical centre ${\mathit{O}}_{\mathcal{C}}$ and which is orthogonal to the camera axis. Consequently, the projection of ${\mathit{O}}_{\mathcal{P}}$ in the image plane is a point at infinity (an ideal point) and the projections of the generatrices all belong to lines parallel to each other.

The 3D curve $\mathcal{C}$ is the intersection between the cone and the model. The set of 3D points on this curve is thus obtained from the intersections of the generatrices of the cone with the model, and they are shown in black in Figure 16. The number of 3D points is equal to the number of generatrices chosen to represent the cone. These 3D points are then projected onto the image to obtain a set of 2D points belonging to the curve $curveProjCone$. These 2D points are shown in black in Figure 16.

The idea now is to obtain the 3D points of the curve $\mathcal{C}$ from the 2D points of the curve ${\mathcal{C}}^{\prime }$. We will assume that the parameters of the cone and the camera are known.

Calculating the Generatrices and the Lines Containing Their Projection

The first step is to calculate the direction vectors ${\mathit{w}}_{\mathit{A}}$ and ${\mathit{w}}_{\mathit{B}}$ of the generatrices $g_A$ and $g_B$. To do this, we need to calculate the angles ${\theta }_A$ and ${\theta }_B$ using the equation 10.

Once the direction vectors have been obtained, we select two 3D points on $g_A$ and $g_B$, project them into the image and compute the parameters of the lines ${\ell }_{g_A}$ and ${\ell }_{g_B}$. The two lines are defined by their respective unit direction vectors ${\mathit{u}}_{\mathit{A}}$ and ${\mathit{u}}_{\mathit{B}}$ and by ${\mathit{A}}_{\mathbf{1}}^{\prime }$ and ${\mathit{B}}_{\mathbf{1}}^{\prime }$, two points belonging respectively to the two lines.

Determining the Points ${\mathit{A}}^{\prime }$ and ${\mathit{B}}^{\prime }$

We wish to obtain the points ${\mathit{A}}^{\prime }$ and ${\mathit{B}}^{\prime }$ which are theoretically the points of tangency of the lines ${\ell }_{g_A}$ and ${\ell }_{g_B}$ with the curve ${\mathcal{C}}^{\prime }$. We will only present the approach for the point ${\mathit{A}}^{\prime }$ since it will be the same for ${\mathit{B}}^{\prime }$.

A first solution would be to say that ${\mathit{A}}^{\prime }$ is the point belonging to the line ${\ell }_{g_A}$ closest to the curve ${\mathcal{C}}^{\prime }$. However, we have chosen a slightly more generic approach, the result of which is illustrated in Figure 17. It consists of considering ${\mathit{A}}^{\prime }$ as the point belonging to the line ${\ell }_{g_A}$ which minimises the sum of the squared distances between itself and the x% of the points on the curve ${\mathcal{C}}^{\prime }$ closest to the line ${\ell }_{g_A}$. Let t be the parameter of the line. We are therefore looking for t such that :

$$t=\underset{t}{argmin}\left\{\sum _{i=1}^n\frac{{({\mathit{A}}^{\prime }-{\mathit{X}}_{\mathit{i}}^{\prime })}^T({\mathit{A}}^{\prime }-{\mathit{X}}_{\mathit{i}}^{\prime })}{h\left({\mathit{X}}_{\mathit{i}}^{\prime }\right)}\right\}\mathrm{with}{\mathit{A}}^{\prime }={\mathit{A}}_{\mathbf{1}}^{\prime }+t{\mathit{u}}_{\mathit{A}}$$

(14)

Where :

• ${\mathit{X}}_{\mathit{i}}^{\prime }$ are the x% of the points on the curve ${\mathcal{C}}^{\prime }$ closest to the line ${\ell }_{g_A}$.
• $h:{\mathbb{R}}^3\to {\mathbb{R}}^{+}$ is a function whose aim is to reduce the impact of the points ${\mathit{X}}_{\mathit{i}}^{\prime }$ far from the line ${\ell }_{g_A}$ in minimisation.

Figure 17 shows the points ${\mathit{A}}^{\prime }$ and ${\mathit{B}}^{\prime }$ obtained, with 40% of the points closest to ${\ell }_{g_A}$ in yellow and 40% of the points closest to ${\ell }_{g_B}$ in orange. Developing the sum to be minimised in 14 gives us :

$$\begin{array}{cc}\hfill \sum _{i=1}^n\frac{{({\mathit{A}}^{\prime }-{\mathit{X}}_{\mathit{i}}^{\prime })}^T({\mathit{A}}^{\prime }-{\mathit{X}}_{\mathit{i}}^{\prime })}{h\left({\mathit{X}}_{\mathit{i}}^{\prime }\right)}& =\sum _{i=1}^n\frac{{({\mathit{A}}_{\mathbf{1}}^{\prime }+t{\mathit{u}}_{\mathit{A}}-{\mathit{X}}_{\mathit{i}}^{\prime })}^T({\mathit{A}}_{\mathbf{1}}^{\prime }+t{\mathit{u}}_{\mathit{A}}-{\mathit{X}}_{\mathit{i}}^{\prime })}{h\left({\mathit{X}}_{\mathit{i}}^{\prime }\right)}\hfill \\ \hfill & =\sum _{i=1}^n\frac{{(t{\mathit{u}}_{\mathit{A}}+{\mathit{v}}_{\mathit{i}})}^T(t{\mathit{u}}_{\mathit{A}}+{\mathit{v}}_{\mathit{i}})}{h\left({\mathit{X}}_{\mathit{i}}^{\prime }\right)}\mathrm{with}{\mathit{v}}_{\mathit{i}}={\mathit{A}}_{\mathbf{1}}^{\prime }-{\mathit{X}}_{\mathit{i}}^{\prime }\hfill \\ \hfill & =\sum _{i=1}^n\frac{(t^2{{\mathit{u}}_{\mathit{A}}}^T{\mathit{u}}_{\mathit{A}}+2t{{\mathit{u}}_{\mathit{A}}}^T{\mathit{v}}_{\mathit{i}}+{{\mathit{v}}_{\mathit{i}}}^T{\mathit{v}}_{\mathit{i}})}{h\left({\mathit{X}}_{\mathit{i}}^{\prime }\right)}\hfill \\ \hfill & =t^2{{\mathit{u}}_{\mathit{A}}}^T{\mathit{u}}_{\mathit{A}}\sum _{i=1}^n\frac{1}{h\left({\mathit{X}}_{\mathit{i}}^{\prime }\right)}+t2\sum _{i=1}^n\frac{{{\mathit{u}}_{\mathit{A}}}^T{\mathit{v}}_{\mathit{i}}}{h\left({\mathit{X}}_{\mathit{i}}^{\prime }\right)}+\sum _{i=1}^n\frac{{{\mathit{v}}_{\mathit{i}}}^T{\mathit{v}}_{\mathit{i}}}{h\left({\mathit{X}}_{\mathit{i}}^{\prime }\right)}\hfill \\ \hfill & =a{t}^2+2bt+c\hfill \end{array}$$

avec

$$\left\{\begin{array}{c}\begin{array}{cc}\hfill & a={{\mathit{u}}_{\mathit{A}}}^T{\mathit{u}}_{\mathit{A}}\sum _{i=1}^n\frac{1}{h\left({\mathit{X}}_{\mathit{i}}^{\prime }\right)}\hfill \\ \hfill & b=\sum _{i=1}^n\frac{{{\mathit{u}}_{\mathit{A}}}^T{\mathit{v}}_{\mathit{i}}}{h\left({\mathit{X}}_{\mathit{i}}^{\prime }\right)}\hfill \\ \hfill & c=\sum _{i=1}^n\frac{{{\mathit{v}}_{\mathit{i}}}^T{\mathit{v}}_{\mathit{i}}}{h\left({\mathit{X}}_{\mathit{i}}^{\prime }\right)}\hfill \end{array}\hfill \end{array}\right.$$

Since a must be positive, the polynomial $a{t}^2+2bt+c$ has a minimum when its derivative cancels, i.e. when :

$$2at+2b=0\Rightarrow t=-\frac{b}{a}$$

This gives :

$${\mathit{A}}^{\prime }={\mathit{A}}_{\mathbf{1}}^{\prime }-\frac{b}{a}{\mathit{u}}_{\mathit{A}}$$

(15)

Using the same procedure, we can obtain the point ${\mathit{B}}^{\prime }$. This produces the line $\left({\mathit{A}}^{\prime }{\mathit{B}}^{\prime }\right)$ (see Figure 17) which will be used to divide the curve ${\mathcal{C}}^{\prime }$ in two.

Separation of the Curve ${\mathcal{C}}^{\prime }$ for Intersection Selection

The points ${\mathit{A}}^{\prime }$ and ${\mathit{B}}^{\prime }$ are determined in a perfectible way in reality. We have therefore defined a minimum distance around the line $\left({\mathit{A}}^{\prime }{\mathit{B}}^{\prime }\right)$ called e which the points of the curve ${\mathcal{C}}^{\prime }$ must respect. Points whose distance orthogonal to the line is less than e are considered indeterminate and will therefore not be taken into account in the 3D reconstruction of the contour. In practice, the value of e will depend on the quality of the measurements taken.

This separation of the ${\mathcal{C}}^{\prime }$ curve into two curves is illustrated in Figure 18 with ${{\mathcal{C}}^{\prime }}_1$ in cyan and ${{\mathcal{C}}^{\prime }}_2$ in magenta.

To obtain ${{\mathcal{C}}^{\prime }}_1$ or ${{\mathcal{C}}^{\prime }}_2$, we had to use the point ${\mathit{p}}_{\mathit{ref}}$ (see equation 12) obtained via the angle $\nu$ (see equation 13) and apply the following two conditions:

$$\left\{\begin{array}{c}\begin{array}{cc}\hfill & {\mathit{X}}^{\prime }\in {{\mathcal{C}}^{\prime }}_2\mathrm{if}sign\left({{\mathit{n}}_{{\mathit{A}}^{\prime }{\mathit{B}}^{\prime }}}^T({\mathit{X}}^{\prime }-{\mathit{A}}^{\prime })\right)=\phantom{-}sign\left({{\mathit{n}}_{{\mathit{A}}^{\prime }{\mathit{B}}^{\prime }}}^T({\mathit{p}}_{\mathit{ref}}-{\mathit{A}}^{\prime })\right)\hfill \\ \hfill & {\mathit{X}}^{\prime }\in {{\mathcal{C}}^{\prime }}_1\mathrm{if}sign\left({{\mathit{n}}_{{\mathit{A}}^{\prime }{\mathit{B}}^{\prime }}}^T({\mathit{X}}^{\prime }-{\mathit{A}}^{\prime })\right)=-sign\left({{\mathit{n}}_{{\mathit{A}}^{\prime }{\mathit{B}}^{\prime }}}^T({\mathit{p}}_{\mathit{ref}}-{\mathit{A}}^{\prime })\right)\hfill \end{array}\hfill \end{array}\right.$$

with ${\mathit{n}}_{{\mathit{A}}^{\prime }{\mathit{B}}^{\prime }}$, the normal vector to the line $\left({\mathit{A}}^{\prime }{\mathit{B}}^{\prime }\right)$.

3D reconstruction

We now know which points ${\mathit{X}}^{\prime }$ on the curve ${\mathcal{C}}^{\prime }$ belong to the curve ${{\mathcal{C}}^{\prime }}_1$ or the curve ${{\mathcal{C}}^{\prime }}_2$. We can therefore apply what was presented in the Section 3.3 and Section 3.5 to reconstruct these points in 3D.

Figure 19 takes the 3D scene from Figure 16 and adds the reconstructed 3D points. The calculated differences between these reconstructed 3D points and the initial 3D points of our simulation (obtained via the intersection between the generatrices of the cone and the gallery model) are below ${10}^{-8}$m, which is very small. These deviations are due to numerical errors in the calculations at the various stages. This proves that these reconstructed 3D points are the same as the initial 3D points and that our method therefore works. We can now test the method with real data.

4. Experiments in a Waterless Aqueduct

The local heritage offers us access, between two of our laboratory buildings, to the Saint-Clément aqueduct (more commonly known as the Arceaux aqueduct) built in the 18th century. It is no longer in use, but in dry periods it is an ideal experimental platform. The gallery is accessible via trapdoors, which once closed plunge you into total darkness. This confined space recreates experimental conditions similar to those in a karstic environment, but in a dry environment whose dimensions are easy to measure. Indeed, the aqueduct is simply a long corridor whose left and right walls are parallel and separated by a width of 62cm, which we measured with a laser rangefinder.

4.1. Our System

The camera used is the Nikon D7000 reflex and we set it to a focal length of 18 mm and a resolution of $2464\times 1632$.

The projector is the DIVEPRO M35 dive light with an aperture angle of ${145}^{\circ }$ in air and ${90}^{\circ }$ in water. As the aperture angle is too large to capture the full projection of light in our experiment, a custom-made plastic cylinder (obtained from the 3D printer) was added to reduce it, as can be seen in Figure 21. Its internal diameter is the same as the diameter of the projector, i.e. $d=$ 53,8mm. The cylinder protrudes from the projector by $l=93$ mm (see Figure 20). We can therefore calculate the half-angle of aperture $\alpha$ as a function of the lengths l and d :

$$\alpha =arctan\left(\frac{d/2}{l}\right)=arctan\left(\frac{26.9}{93}\right)={16.13}^{\circ }$$

We will use this value to check the order of magnitude of future $\alpha$ values when we estimate the cone.

The camera and projector are placed on a rigid support as shown in Figure 21. They had to be fixed to the support because the slightest rotation of one in relation to the other alters the transformation between the camera and the projector, which is estimated during the projector calibration.

To get an idea of the magnitude of this transformation, we used a caliper to measure the translations between the camera and the projector along the camera’s x, y and z axes. The rotation between the two is supposed to be small because we have chosen an identical orientation for both. Here is a summary of the various measurements of the projector parameters (lengths in m):

$$\left\{\begin{array}{c}\begin{array}{cc}\hfill \alpha & \approx {16.13}^{\circ }\hfill \\ \hfill {{\mathit{O}}_{\mathcal{P}}}_{\left\{\mathcal{C}\right\}}& \approx {\left[\begin{array}{ccc}0.2& 0& -0.07\end{array}\right]}^T\hfill \\ \hfill ∥{{\mathit{O}}_{\mathcal{P}}}_{\left\{\mathcal{C}\right\}}∥& \approx 0.212m\hfill \\ \hfill {\mathit{d}}_{\left\{\mathcal{C}\right\}}& \approx {\left[\begin{array}{ccc}0& 0& 1\end{array}\right]}^T\hfill \end{array}\hfill \end{array}\right.$$

4.2. Camera Calibration Results

To calibrate our camera, we use the method of [30] with the help of eight images of a flat chessboard test pattern measuring 7x10 squares and 36mm wide (see Figure 22). In each of the eight images, the position of the test pattern relative to the camera is different, as can be seen in Figure 23, which is a 3D representation of the test pattern in eight different positions. Here are the camera parameters estimated via calibration:

• the intrinsic matrix K obtained is :

  $$K=\left[\begin{array}{ccc}1908.56& 0.0& 1227.53\\ 0.0& 1909.94& 832.48\\ 0.0& 0.0& 1.0\end{array}\right]$$
• the resulting focal length is 18.2 mm,

It can be seen that the orders of magnitude of the various camera parameters are consistent with the manufacturer’s parameters. In addition, the average reprojection errors are of the order of 1/10 of a pixel. This level of error is quite satisfactory for all eight shots, each containing 56 calibration points.

4.3. Projector Calibration Results

To calibrate the projector, i.e. estimate the parameters of the cone, we used a white wall on which a chessboard pattern was hung. The first step is to obtain 3D points belonging to the cone by detecting the contour of the light projected onto the wall, as explained in the section entitled 3.2. To maximise the contrast of the light projected onto the white wall, the only light source in the room is the projector.

For this calibration, we captured five images of light projected onto the wall with the axis of the cone almost orthogonal to the wall (see Figure 24). The shape of the projected light therefore resembles a circle.

Each detected contour is transformed into an ellipse, shown in blue in Figure 24. To obtain the ellipses we use the method detailed in [28] but iteratively with the RANSAC algorithm [31] to eliminate certain outlying contour points. The contour points eliminated by RANSAC are shown in red in Figure 24 while the points used to obtain the ellipses are in green.

Using the pattern detected in each image, we obtain the equations for the five planes. From these five planes and the five ellipses obtained, the five elliptical sections are calculated. They are represented in 3D in Figure 25 with the camera and the test pattern.

The image on the left shows all the elements expressed in the chessboard pattern frame, so you can see the different camera positions used to capture the five images.

The image on the right shows all the elements expressed in the camera frame. It is in this frame of reference that the elliptical sections must be expressed in order to estimate the cone.

We extract 3D points expressed in the camera frame of reference for each elliptical section. In our case, we recovered 200 points per section, making a total of 1000 3D points. Note that these 3D points are at a distance from the camera frame of between 1,50m and 2,17m.

Now that we have obtained 1000 3D points that are supposed to belong to the cone, the next step is to define an initial cone for our iterative minimisation algorithm. This is shown in red in Figure 26. We chose an initial cone that is far from the solution for greater visibility in the figure and also to prove that the method is capable of converging even when the initial cone is far from the solution. Obviously, this is an experimental result and not a mathematical proof. This optimisation problem is probably non-convex. Ideally, therefore, an initial cone close to the solution should be chosen using the aperture angle previously calculated and the relative pose measurement between the projector and the cone.

We therefore apply our algorithm and obtain the optimised cone shown in grey. Its estimated parameters (lengths in metres) are :

• $\alpha ={14.79}^{\circ }$
• ${{\mathit{O}}_{\mathcal{P}}}_{\left\{\mathcal{C}\right\}}={\left[\begin{array}{ccc}0.1956& 0.0092& -0.0759\end{array}\right]}^T$
• $∥{{\mathit{O}}_{\mathcal{P}}}_{\left\{\mathcal{C}\right\}}∥=0.210m$
• ${\mathit{d}}_{\left\{\mathcal{C}\right\}}={\left[\begin{array}{ccc}-0.013& 0.080& 0.997\end{array}\right]}^T$

To ensure the validity of the cone obtained, we begin by comparing the estimated parameters with the previously measured parameters. The Table 1 compares the estimated opening angle and the estimated gap with the measurements. As our measurements can be improved, this comparison at least proves that the order of magnitude of the estimated parameters is consistent.

Another point to check is the minimisation error for estimating the cone. In our case, this error is directly linked to the orthogonal distances between the 3D points of the elliptical sections and the cone, since it is the root mean square (RMS) of these distances that we are trying to minimise during estimation. In Figure 26, the colour of each 3D point is defined by its orthogonal distance from the cone via the colormap Jet (the colormap bar is shown in the figure) and the following function :

$$\begin{array}{cc}\hfill f:{\mathbb{R}}^{+}& \to [0,1]\hfill \\ \hfill d& \mapsto \frac{d}{d+0.1}\hfill \end{array}$$

with d the orthogonal distance. m is chosen as the median of all these distances (here $m=1.4$mm).

The statistics for these distances can be found in the Table 2. In particular, we have a maximum distance of 4,2mm and an overall mean of 1,5mm. The overall RMS is 1,9mm which corresponds to approximately 0,13% of 1,50m which is the smallest distance between a 3D point of an elliptical section and the camera frame.

In view of the analysis of the orthogonal distances, the minimisation is satisfactory, confirming once again that the estimated parameters of the cone are at least consistent with their true values.

4.4. 3D Results

The shape of our aqueduct is a long narrow corridor, which makes it possible to fully distinguish the light projection in our images. We captured six images taken at different distances from a chessboard pattern placed in the aqueduct. For each image, the light contours were checked manually, as we wanted to isolate the triangulation method from the contour detection method, which is still in need of improvement. For each contour extracted, we manually delimited part of the contour points on the left wall and another part on the right wall. These delimitations can be seen in Figure 27. It shows an example of an extracted contour in two images, with contour points on the left wall in red, on the right wall in blue and the rest of the contour points in green.

Only the contours of the left and right walls will be used for triangulation, as this will allow us to check that the 3D points obtained respect the dimensions of the aqueduct corridor. To obtain these 3D points, it is necessary, as usual, to define which intersection is relevant for each radius associated with the contour points. In our situation, we can see that each contour point on the left wall implies a first intersection and each contour point on the right wall a second intersection. Figure 28 illustrates this for the two images in Figure 27.

Figure 29 shows the result of the reconstruction of the contours of the six images. The shape of the 3D reconstructions of the contours is close to a hyperbola, which is consistent since the intersection between a plane and a cone in this configuration is supposed to be a hyperbola.

We’re now going to compare this 3D reconstruction with reality. The aqueduct is a corridor whose left and right walls are parallel and separated by a width of 62cm. We can therefore check that :

• the 3D points resulting from the left part (respectively right part) of the contour belong to the same plane so they are supposed to be coplanar,
• If we estimate a plane from the left 3D points and a plane from the right 3D points, they must be parallel and separated by a distance of approximately 62cm.

As shown in Figure 29b, we get estimations of the left wall in red and the right wall in blue by using a least square method. To evaluate coplanarity, we compute the deviation between 3D points and their associated plane. For left points, we get an average deviation of 13,9mm and for right points 24,1mm (see Table 4).

To check the gap between the two estimated planes, we compute the distances between the left 3D points (respectively the right 3D points) and the right plane (respectively the left plane). The average distance between the left points and the right plane is 584,9mm. The average distance between the right points and the left plane is 583,6mm.

The two estimated planes are almost parallel since the angle calculated between the 2 normal vectors is 178,9°.

5. Discussion

Authors should discuss the results and how they can be interpreted from the perspective of previous studies and of the working hypotheses. The findings and their implications should be discussed in the broadest context possible. Future research directions may also be highlighted.

6. Conclusions

The aim of this paper was to propose a 3D reconstruction solution designed for underwater galleries found in karstic environments. The study of these environments represents a considerable challenge for the future, as they could provide part of our water needs. Exploration remains the best approach for acquiring reliable data on the structure, and is mainly carried out by qualified divers using topography methods with Ariane wires and section-by-section surveys. But on extensive karstic networks, this approach presents a major risk, as the duration of these surveys is so time-consuming, leading divers to long stops. The future is therefore more likely to lie with robotics, but there are still many technical challenges to be overcome before underwater drones that are sufficiently autonomous can be sent out to explore this type of environment. Our aim was to propose an original method of 3D scanning in this type of environment. Our approach uses a combination of a camera and a cone-shaped projector. The parameters of the light cone have been fixed: its centre, its direction vector and its angular aperture. A method for calibrating the cone was presented, enabling the position of the 3D points obtained by intersecting the visual rays of the contour of the projector halo as seen by the camera with the cone whose parameters have been estimated. The calibration of the camera is done thanks to a classical approach. A simulation of this method was presented and results validate the proposal. A first experimentation in real conditions was conducted in a disused aqueduct (without water) which had the advantage of reproducing an environment close to that of a karstic aquifer: narrow, without light and above all without water to continue our evaluation. The camera calibration is validated with an error reconstruction less than 0.1 pixel and camera parameters close to technical values (see 4.2). The cone calibration leads to an angles error less than 2 degrees and distance with respect to the camera less than 2 mm. The results were validated by testing the coplanarity of the aqueduct walls to within 25mm and estimating the known distance between these walls to within 36mm (see Table 4). Experiments have be done in real karstic environment but due to complex meteorological conditions, visual acquisition was not possible. Thus, we aim to begin in a more controlled environment : a swimming pool and then finish in a real aquifer with clearer water.

To explore a new configuration of this system, we began to make some simulations with one projector and four cameras (Figure 30) to solve the error of the transition zone (see Figure 19). We also think about using several projectors with localisation sensors with respect to one or more cameras to have different configurations and measurement redundancy.