L Luneville
email: laurence.luneville@cea.fr

D Simeone
email: david.simeone@cea.fr

Impact of the ageing of materials on primary damage induced by different neutron spectra

Keywords: ageing, neutron-isotope cross section, primary and recoil spectra

The recent nuclear evaluations like the Thalys based Evaluation Nuclear Data Library (TENDL) are now able to accurately estimate either from Monte Carlo simulations or experimental results the many elastic and inelastic neutronatoms interactions and allow calculating more realistically primary damage induced by nuclear reactions. Even if these calculations do not take into account relaxation processes occurring at the end of the displacement cascade, they can accurately describe primary and recoil spectra in different reactors or ion beam accelerators. Since neutrons are sensitive to isotopes and not to chemical elements, these spectra must be calculated taking into account all neutron-isotopes interactions separately and then be summed over the relative abundance of different isotopes forming chemical elements. Such a procedure is applied in the DART code developed to compare irradiation in different facilities . Such a procedure gives an accurate estimation of primary damage induced by different fluxes within the Binary Collision Approximation using the TENDL-2015 library allowing to mimic primary damage produced in different nuclear facilities with ion beams. In this work we discuss the potentiality of such an approach computing primary damage in Ni based steel springs of nuclear reactors for different irradiation times, pointing out the importance of ageing, i.e the modification of isotopic abundance of Ni isotopes with the irradiation fluence, on these spectra. This calculation taking into account accurate neutron isotope cross sections is crucial to determine the nature and the energy of ions used in ion beam facilities to mimic the micro structural evolution of actual and future materials for nuclear plants.

Introduction

When a material is subjected to a neutron flux, atoms of the solid can be knocked from their equilibrium positions with a wide range of energies. The nature and the energy of recoil atoms produced during inelastic neutron isotopes collisions determine for a large part the primary damage responsible for the long term stability of the material under irradiation. Averaging these recoils energies over the isotopic composition allows studying the impact of ageing, i.e. the variation of the isotopic composition versus the fluence, on the primary damage. The detailed description of calculation performed by the DART code 1 has already been discussed in detail in previous articles [START_REF] Luneville | Calculation of radiation damage induced by neutrons in compound materials[END_REF] and only specific applications of this code to study the ageing of materials are discussed in this work. To quantify the effect of ageing of materials, different Primary Knocked Atoms (PKA) spectra and displacement cross sections [START_REF] Odette | Neutron-Energy-Dependent Defect Production Cross Sections for Fission and Fusion Applications[END_REF][START_REF] Simeone | Contribution of recoil atoms to irradiation damage in absorber materials[END_REF] have been computed in nickel-containing alloys with different compositions for different facilities.

The impact of the isotopic composition on the primary damage on Ni element is discussed in detail in the first part of this work pointing out the need to have accurate neutron isotopes cross section describing all possible open channels. This last point implies an large and accurate estimation of all stable and non stable nuclei produced during a neutron isotope interaction. All calculations presented in this work were performed using the TENDL-2015 library which combines Monte Carlo and analyses of experimental measurements to compute these data [START_REF] Koning | TALYS-1.0[END_REF]. A close inspection of all neutron iso-tope cross sections points out that primary damage, i.e. displacement of atoms induced by the neutron impinging over a time scale smaller than 10 -9 second, is only produced by 58 Ni and 59 Ni isotopes. The dependence of the PKA spectra and the displacement cross sections as a function of the different isotopes and the neutron fluxes of different facilities is discussed in the second part of this work, allowing then to discuss the effect of aging in Ni based steels.

Isotopic composition and primary damage

When a neutron interacts with an atom at rest in a crystal, a part of the kinetic energy of the neutron is given to the atom leading to the recoil of the atom (elastic cross section) or a nuclear reaction (inelastic cross section). The cross sections associated with a specific interaction determine the intensity of the process and the part of the energy given to the atom at rest (elastic collision) or to daughter atoms produced during the nuclear reaction (inelastic collision). Obviously, the cross sections associated to each open channel is different for each isotope of a given chemical element. To accurately compute the all possible nuclear interactions produced by impinging neutrons in a given material, all neutron isotopes cross sections must be computed and then averaged according to the isotopic composition of each chemical element forming the crystal. This procedure is applied in the DART code [START_REF] Luneville | Calculation of radiation damage induced by neutrons in compound materials[END_REF] and such a procedure is different from those used in other codes like SPECTER [START_REF] Greenwood | SPECTER: neutron damage calculation for materials irradiations[END_REF] which computes only neutron-element cross sections. Such a procedure thus allows following the evolution of primary damage versus the material ageing, i.e. the variation of its isotopic composition as a function of the dose.

To point out the impact of the isotopic composition on the primary damage for different facilities, we compute primary spectra associated with the energy redistribution of primary knocked on atoms for different facilities and the displacement cross sections versus the energy of the incident neutron, at different isotopic values for pure nickel atoms. The variation of the isotopic composition of 58 Ni versus the irradiation time may give some clues to explain the unexpected embrittlement of Ni based alloys springs observed in Candu reactors [START_REF] Griffiths | The effect of irradiation on Ni-containing components in Candu reactor cores: a review[END_REF][START_REF] Rowcliffe | Perspectives on radiation effects in nickel-base alloys for applications in advanced reactors[END_REF] Figure 1 displays the neutron spectra of different nuclear facilities used to study the aging of Ni based alloys. To highlight the effect of inelastic collisions, PKA spectra and displacement cross sections induced by different ion beams (100 keV Xe and 800 keV Ar) were also calculated for pure Nickel.

Nuclear cross sections of Ni isotopes

The pure nickel is mainly composed by 3 different isotopes. The 58 Ni is the most abundant isotope (67.8%) followed by 60 Ni (26.2%) and finally the 62 Ni (3.7%). The two last isotopes 60 Ni and 62 Ni are stable and their inelastic neutron isotopes cross section can be neglected. Even if neutron-isotope cross sections are similar for all three isotopes, the 58 Ni is the most abundant isotope and is the only one able to produce the 59 Ni metastable nucleus according to the capture reaction:

58 Ni + n → 59 Ni (1)
The 58 Ni capture cross section follows a square root law of the neutron kinetic energy in the center of mass frame. This capture cross section varies from 100 barns to 0.01 barns for a neutron kinetic energy ranging from 10 -5 eV to 10 3 eV.

The unusual properties of metastable 59 Ni isotope is well established [START_REF] Griffiths | The effect of irradiation on Ni-containing components in Candu reactor cores: a review[END_REF][START_REF] Greenwood | A new calculation of thermal neutron damage and helium production in nickel[END_REF]:

59 Ni + n → 56 Fe + α Q = 5.098MeV (2) 59 Ni + n → 59 Co + H Q = 1.85MeV (3) 59 Ni + n → 60 Ni + γ Q = 11.38MeV (4)
The 59 Ni(n, α) reaction releases 5.1 MeV leading to the emission of a 4.8 MeV α particle and 340 keV Fe recoils for thermal neutrons. The second reaction 59 Ni(n, p) mainly generates a 1.85 MeV proton and 31 keV 59 C recoils for thermal neutrons. The last nuclear reaction only induces the emission of a γ which does not participate to the production of primary damage.

Energy spectra of recoil atoms

Neglecting the γ emission, the number of daughter recoils produced during the nuclear reaction depends on the n-59 Ni capture cross section versus the kinetic energy of the impinging neutrons. Such capture cross section for two distinct channels are given by neutron libraries. In this work, all n-59 Ni capture cross section were computed from the TENDL-2015 [START_REF] Koning | TALYS-1.0[END_REF] cross section data library which contains the most accurate evaluations on 59 Ni. Obviously, the energy of recoil atoms ( 56 Fe and 59 Co in our case) is not only a function of the kinetic energy of the incident neutron but also of the angular distribution of emitted particle in the center of mass frame:

T (µ, E) = 1 (A + 1) (E * -2µ A E * E a + A E a ) (5)
where

E * = A+1-A 1+A E and E a = Q + A A+1 E
in the center of mass frame. In these expressions, A is the mass ratio of incident particles and A' is the mass ratio of ejected particles [START_REF] Macfarlane | The NJOY nuclear data processing system[END_REF]. Then it appears that the energy of recoils induced by nuclear reaction listed in Eq. ( 2) depend on the kinetic energy of incident neutron and then the neutron fluxes of different facilities.

This expression points out that the energy of recoil atoms are function of the kinetic energy of the incident neutron and then PKA spectra induced by different facilities will be function of their neutron flux plotted on Fig. 1. To point out the impact of neutron fluxes on PKA spectra, different facilities associated with fusion aznd fusion reactors are plotted on this graph. In order to achieve the comparison with ion beam facilities, typical ion beam spectra which behave as Dirac distributions are also plotted on the same graph just for comparison.

PKA spectra and displacement cross sections

To highlight the impact of different Ni isotopes on the PKA spectra, Figure 2 displays the variation of the PKA spectra as a function of its kinetic energy for 58 Ni and 59 Ni isotopes irradiated in the HFIR reactor [START_REF] Blosser | Neutron Flux and Neutron and Gamma-Ray Spectra Measurements at the HFIR[END_REF] used to test materials for fission reactors,.e.i exhibiting a large thermal neutron flux (cf Fig. 1). Two distinct peaks clearly appear on the PKA distribution for the 59 Ni isotope. 59 Ni PKA distribution pointing out the impact of the Ni isotopes on the primary damage production. These calculations result from the integration of neutronisotopes cross sections with the neutron spectrum of the HFIR facility. All these calculations were performed using the JENDL-4.0 library.

Detailed study of nuclear reaction in the 59

Ni isotope In order to investigate the origin of the two peaks observed on the 59 Ni PKA spectra plotted on Figure 2, PKA spectra associated with the elastic collision, the (n, α) and the (n,p) nuclear reactions were plotted on Figure 3. The analysis of these graphs clearly points out that peaks observed on the 59 Ni PKA spectrum are associated with the recoil of 56 Fe and 59 Co and their intensities and positions are function of the incident neutron flux.

In order to estimate the impact of the recoil on the formation of primary damage production, the displacement cross section computed within the BCA framework was computed for different nuclear reactions for the HFIR facility in 59 Ni isotope. Figure 4 displays the variation of the displacement cross section for each nuclear reaction as a function of the kinetic energy of the incident neutron in the laboratory frame.

The analysis of Figure 4 clearly points out that displacement cross section induced by the 56 Fe and 59 Co recoil are dominant over a large energy range. These displacement cross sections are one order of magnitude larger than displacement cross section for the 58 Ni isotope. From the computation of the partial displacement cross section we can conclude that the main part of primary damage produced in Ni based materials results from the recoil atoms induced by nuclear reaction on the 59 Ni isotope. This analysis illustrates the need to compute isotope by isotope the primary damage in materials in order to accurately account all possible nuclear reactions.

Aging of Ni based materials for different facilities

Despite the displacement cross section for 59 Ni isotope is one order of magnitude larger than the one for 58 Ni isotope as pointed on Figure 4, the primary damage production rate (dpa/s) calculated within the BCA framework results from the integration of these displacements over the different neutron fluxes for a given nuclear facility.

To study in detail the effect of aging of material submitted to long time neutron flux, primary damage production rates were computed on 10% Ni based Fe steel for different facilities. To mimic the impact of aging on these damage rates, the primary damage production rates were computed for a pristine sample made of 90% of iron and 10% of 58 Ni isotopes and another sample made of 90% of iron and 10% of 59 Ni assuming that all 58 Ni isotopes are transmuted in 59 Ni.

Figure 5 displays the ratio between the primary damage production rates for different facilities. The analysis of this figure points out that the ratio of primary damage production rates is highly sensitive to the thermal neutron fluxes in agreement with the narrow peaks observed on the displacement cross section of the 59 Ni isotope (cf Figure 4).

Despite the appearance of the 59 Ni isotope is responsible for a large increasing of the primary damage pro- duction rate, this rate is mainly controlled by the fraction of thermal neutrons. For HFIR and PWR facilities exhibiting important thermal neutron fluxes (cf Figure 1), this ratio is important. For 800 keV Ar and 100 keV Xe irradiation, this ratio does not evolve because these irradiation do not induce any nuclear reaction in the material. Fusion facilities like ITER and IFMIF exhibiting small thermal neutron fluxes as displayed in Figure 1 are insensitive to the isotopic composition of the steel.

Conclusion

To discuss the importance of the isotopic composition on the primary damage, we study in detail the impact of different neutron isotope cross section in Ni elements, composed of different isotopes. From a careful analysis of inelastic neutron isotope cross section using the TENDL-2015 library, it appears that the main part of PKA and displacement cross sections results from the transmutation of the 59 Ni isotope. The analysis of the displacement cross section, calculated by the DART code which computes PKA and displacement cross section for each isotope, exhibit sharp peaks due to the recoils of 56 Fe and 59 Co atoms resulting from the transmutation of the 59 Ni isotope. Calculations performed with the TENDL-2015 library clearly show that these sharp peaks are centered in the epi-thermal domain, i.e. neutron with kinetic energies varying from few hundredth eV to few eV. From this analysis, it was possible The example studied in this work clearly demonstrates the need to compute primary damage production from the knowledge of the distribution of isotopes in the material and from very accurate neutron-isotope cross section data libraries.

Figure 1 :

 1 Figure 1: Neutron and ions spectra φ for different nuclear reactors and ion beam facilities (pink line: ITER (fusion reactor), turquoise line: HFIR (High Flux Isotope Reactor), navy line: IFMIF (International Fusion Materials Irradiation Facility), green line: PWR (Power Water Reactor); red line: Ar 800 keV; blue line: Xe 100 keV)

Figure 2 :

 2 Figure 2: PKA spectra as a function of the energy of recoil atoms T in pure 58 Ni (magenta line) and 59 Ni (black line) isotopes. Two distinct peaks are clearly observed on the59 Ni PKA distribution pointing out the impact of the Ni isotopes on the primary damage production. These calculations result from the integration of neutronisotopes cross sections with the neutron spectrum of the HFIR facility. All these calculations were performed using the JENDL-4.0 library.

Figure 3 :

 3 Figure 3: Comparison of PKA spectra as a function of the recoil energy T for different nuclear reactions occurring in the impinging of the 59 Ni isotope with a neutron (green line : (n, α) nuclear reaction; blue line: (n,p) nuclear reaction; red line: elastic collision; black line: all processes are taken into account). All these calculations were performed for the HFIR reactor with the TENDL-2015 library on pure 59 Ni isotope.

Figure 4 :

 4 Figure 4: Comparison of the displacement cross sections computed within the BCA framework as a function of the kinetic energy of the incident neutron in the laboratory frame for all nuclear reactions occurring in pure 59 Ni isotope (green line : (n, α) nuclear reaction; blue line: (n,p) nuclear reaction; red line: elastic collision; black line: all processes are taken into account). All these calculation were performed for the HFIR reactor with the TENDl-2015 library.

Figure 5 :

 5 Figure 5: Comparison of the ratio of the primary damage production rate in 10% 59 Ni and 10% 58 Ni based Fe steel (10%Ni and 90% Fe in atomic fraction) for different nuclear facilities. Depending on the neutron spectra of each facility, this ratio can differ from one order of magnitude (red line: ratio equal to unity).

Acknowledgement

We acknowledge Dr. Greenwood for providing us the characteristic neutron fluxes of different nuclear facilities, Dr J.-C. Sublet for giving us the TENDL-2015 library and for his helpful advices on the underlying physics beyond the description of inelastic cross sections neutron-isotopes and all members of the international DPA CRP from IAEA for helpful discussions on primary damage calculations.


Robin Bruneau

Baptiste Brument

Lilian Calvet

Matthew Cassidy

Jean Mélou

Yvain Quéau

Jean-Denis Durou

François Lauze

Multi-view Stereo of an Object Immersed in a Refractive Medium

In this article we show how to extend the multi-view stereo technique when the object to be reconstructed is inside a transparent -but refractive -medium, which causes distortions in the images. We provide a theoretical formulation of the problem accounting for a general, non-planar shape of the refractive interface, and then a discrete solving method. We also present a pipeline to recover precisely the geometry of the refractive interface, considered as a convex polyhedral object. It is based on the extraction of visible polyhedron vertices from silhouette images and matching across a sequence of images acquired under circular camera motion. These contributions are validated by tests on synthetic and real data.

Introduction

Many Natural History Museums have collections of fragile specimens immersed in a transparent medium, for instance a prehistoric insect trapped in a block of amber, or a small animal conserved in a solution of formaldehyde or alcohol in a glass jar (see Figure 1), which constitute an invaluable source of information for evolutionary scientists.

Exploitation of these resources necessitates "see through" 3D estimation techniques. Computerised tomography (CT) provides a 3D reconstruction, however it requires special CT scanners that can handle objects of various sizes. This is expensive and time consuming, and may be feasible for a limited amount of samples but out of reach for large collections (thousands of amber pieces, or linear kilometers of specimens in jars). Photogrammetric 3D scanning is a reasonable alternative, which however requires overcoming several difficulties. First of all, as described by Snell's law, if the refractive medium deflects the trajectory of light rays that pass through it, the shape of the interface, i.e. the separation surface between air and this medium, has a strong influence on the appearance of the object to be reconstructed. In the example on the left in Figure 1, it is not forbidden to cut the amber block to give it a particular shape, as long as this does not alter the insect. The shape of the refractive interface will thus be assumed to be known. However, unlike most existing techniques, we do not restrict our attention to the planar case.

Refraction can cause other difficulties. Because it is wavelength dependent, it causes light dispersion. In addition, an inhomogeneous transparent medium may have a spatially-varying refractive index, resulting in light rays that are not straight. Taking a varying refractive index into account does not present an intractable difficulty, but estimating this index at each point in the medium is in itself a difficult problem. We will therefore assume that the refractive medium is homogeneous and that the refractive index varies little with wavelength.

Contributions. We introduce a general 3D reconstruction method for objects trapped inside a homogeneous refractive medium with arbitrary, known interface, from a set of calibrated multi-view stereo images. The paper is organized as follows. In Section 2, we review the main studies taking refraction into account. We show in Section 3 how to extend the multi-view stereo technique in the presence of an interface. The main difficulty of this extension consists in predicting the projection in an image of a 3D point located in the refractive medium. A first set of tests on synthetic images is conducted in Section 4 to validate this method. In Section 5, we show how to validate it on real data, which requires additional information about the 3D scene. Finally, we conclude in Section 6 by mentioning some possible extensions. The end-to-end pipeline described in this paper extends two of our earlier papers [START_REF] Cassidy | Refractive Multi-view Stereo[END_REF][START_REF] Brument | A Shape-from-silhouette Method for 3D-reconstruction of a Convex Polyhedron[END_REF].

Related work

Refraction in computer vision has been considered in various ways: it can be considered as a bias to correct in order to use classic vision techniques, purposely used in an active system, or explicitly modeled in a refractive 3D reconstruction pipeline.

Refraction compensation in classic vision. Light rays from a point source naturally tend to diverge. By refraction, a lens allows to deviate their trajectories to make them converge in an image point. Refraction is therefore the key to optical instruments, with the caveat that manufacturing a photographic lens assembly requires very precise alignment of the lenses that make it up, in order to limit undesired effects, called aberrations. However, if the scene itself contains transparent, i.e. refractive, objects, this modifies the appearance of opaque objects, which appear distorted. A number of papers have addressed the correction of these distortions when the transparent object, such as a window pane, is attached to the front of the camera, in which case the distortions can be calibrated, and standard 3D reconstruction pipelines can then be employed. For example, while it is relevant to use photogrammetry to study convection in a water-filled tank, Maas shows in [START_REF] Maas | New developments in multimedia photogrammetry[END_REF] how to take into account refraction through both sides of the glass to improve measurement accuracy. In [START_REF] Łuczyński | Image rectification with the pinax camera model in underwater stereo systems with verged cameras[END_REF], Łuczyński et al. show how to correct the images acquired by a pair of calibrated underwater cameras, protected by a glass plate, to restore epipolar geometry. The same idea of pre-correcting the images has also been explored in [START_REF] Agrafiotis | Correcting image refraction: Towards accurate aerial image-based bathymetry mapping in shallow waters[END_REF][START_REF] Wu | Accurate binocular stereo underwater measurement method[END_REF], as well as in [START_REF] Agrafiotis | DepthLearn: learning to correct the refraction on point clouds derived from aerial imagery for accurate dense shallow water bathymetry based on SVMs-fusion with LiDAR point clouds[END_REF] where correction is carried out by neural networks. It has also been shown in [START_REF] Ichimaru | Underwater Stereo Using Refraction-Free Image Synthesized From Light Field Camera[END_REF][START_REF] Zhang | On-site calibration of underwater stereo vision based on light field[END_REF] that light field cameras can provide useful clues for correcting refraction in the images.

Active refraction techniques. On the other hand, several studies, which can be grouped under the term active refraction, take advantage of refraction to allow 3D reconstruction from a single view, thanks to image duplication, using either a bi-prism [START_REF] Lee | A biprism-stereo camera system[END_REF][START_REF] Yamashita | Monocular Underwater Stereo -3D Measurement Using Difference of Appearance Depending on Optical Paths[END_REF], or a glass plate rotating around the optical axis [START_REF] Chen | Depth from refraction using a transparent medium with unknown pose and refractive index[END_REF][START_REF] Gao | A Refractive Camera for Acquiring Stereo and Super-resolution Images[END_REF].

Estimation of a refractive interface. In many cases, the surface of a transparent object (the refractive interface) may be needed. In the case of water, Morris uses a pattern refracted through water surface in [START_REF] Morris | Image-based water surface reconstruction with refractive stereo[END_REF], Morris and Kutulakos propose to map points of interest seen through transparency [START_REF] Morris | Dynamic Refraction Stereo[END_REF]. Another approach is that of Ben-Ezra and Nayar [START_REF] Ben-Ezra | What does motion reveal about transparency?[END_REF], who fit the parameters of a surface model in order to reproduce the distorted image of an object of known geometry. The 3D reconstruction of transparent objects has also been revisited recently under the neural network (NN) perspective [START_REF] Li | Through the Looking Glass: Neural 3D Reconstruction of Transparent Shapes[END_REF][START_REF] Shao | Polarimetric Inverse Rendering for Transparent Shapes Reconstruction[END_REF].

Bathymetry. An obvious need for understanding refraction is remote-sensing bathymetry. Its goal is to correct depth measurements in oceans and lake or river beds, which differs somewhat from our reconstruction problem, though uses Snell's law. This has been explored in numerous works, for instance in the work of Murase et al. [START_REF] Murase | A Photogrammetric Correction Procedure for Light Refraction Effects at a Two-Medium Boundary[END_REF], the already mentioned papers by Agrofiotis et al. [START_REF] Agrafiotis | DepthLearn: learning to correct the refraction on point clouds derived from aerial imagery for accurate dense shallow water bathymetry based on SVMs-fusion with LiDAR point clouds[END_REF][START_REF] Agrafiotis | Correcting image refraction: Towards accurate aerial image-based bathymetry mapping in shallow waters[END_REF] but also Woodget et al. in [20] or Cao et al. in [START_REF] Cao | Universal algorithm for water depth refraction correction in through-water stereo remote sensing[END_REF] to cite a few.

Classical framework with refraction adaptation. Classical components of 3D vision systems have been adapted to the presence of a refractive interface: calibration, absolute and relative camera pose estimation and refractive structure-from-motion (RSfM), refractive multi-view stereo (RMVS), refractive photometric stereo (RPS), and refractive multi-view photometric stereo (RMVPS).

Sturm discusses generic camera models for SfM in [START_REF] Sturm | Multi-view geometry for general camera models[END_REF], including the refractive axial camera. Chari and Sturm extends epipolar geometry in [START_REF] Chari | Multi-View Geometry of the Refractive Plane[END_REF] in the case of a camera separated from the 3D scene by a planar interface, using a 12 × 12 fundamental matrix.

A pinhole/axial camera model is proposed by Łuczyński et al. in [START_REF] Łuczyński | Image rectification with the pinax camera model in underwater stereo systems with verged cameras[END_REF] and calibration for it. A fringe projection system is studied by Chen et al. in [START_REF] Chen | Three-dimensional reconstruction from a fringe projection system through a planar transparent medium[END_REF]. Jordt an coworkers study the case of an underwater camera and its consequences in a series of papers. In [START_REF] Jordt-Sedlazeck | Refractive calibration of underwater cameras[END_REF], geometric and radiometric calibration procedures are proposed, and pose by structure-from-motion in [START_REF] Jordt | Refractive 3D reconstruction on underwater images[END_REF][START_REF] Jordt-Sedlazeck | Refractive Plane Sweep for Underwater Images[END_REF], as do Kang et al. [START_REF] Kang | Two-view underwater 3D reconstruction for cameras with unknown poses under flat refractive interfaces[END_REF], Qiao et al. [START_REF] Qiao | Structure from Motion of Underwater Scenes Considering Image Degradation and Refraction[END_REF], Suresh et al. [START_REF] Suresh | Through-water stereo slam with refraction correction for AUV localization[END_REF], Chadebecq et al. [START_REF] Chadebecq | Refractive Two-View Reconstruction for Underwater 3D Vision[END_REF], Hu et al. [START_REF] Hua | Refractive Pose Refinement[END_REF]. Specific pose optimization problems under flat refractive interfaces are discussed in [START_REF] Haner | Absolute pose for cameras under flat refractive interfaces[END_REF]. In most of these papers, validation is carried out on underwater images -a relatively recent survey of underwater 3D imaging techniques can be found in [START_REF] Castillón | State of the Art of Underwater Active Optical 3D Scanners[END_REF]. Complete reconstruction pipelines have been proposed, especially for underwater camera [START_REF] Jordt | Refractive 3D reconstruction on underwater images[END_REF][START_REF] Chadebecq | Refractive Two-View Reconstruction for Underwater 3D Vision[END_REF]. Though the majority of refractive interfaces are planar, series of other situations have also been explored. The case of a camera inside an underwater dome is studied in [START_REF] She | Refractive geometry for underwater domes[END_REF]. More complex situations have explored when seeing aerial objects from underwater through a randomly moving surface in [START_REF] Alterman | 3D in natural random refractive distortions[END_REF][START_REF] Alterman | Triangulation in Random Refractive Distortions[END_REF]. The converse situation (air to water) is explored in [START_REF] Xiong | In-the-wild single camera 3D reconstruction through moving water surfaces[END_REF] using a differentiable/NN framework. Many authors have studied underwater extension of photometric stereo (PS), often focusing on scattering media properties while discarding photometric stereo, or attempting to incorporate via camera calibration as in [START_REF] Tsiotsios | Backscatter compensated photometric stereo with 3 sources[END_REF][START_REF] Murez | Photometric Stereo in a Scattering Medium[END_REF]. Narasimhan et al. incorporate explicitly refraction in PS, among other techniques, in [START_REF] Narasimhan | Structured light in scattering media[END_REF]. Fan et al. [START_REF] Fan | Refractive laser triangulation and photometric stereo in underwater environment[END_REF][START_REF] Fan | Underwater optical 3-D reconstruction of photometric stereo considering light refraction and attenuation[END_REF] use it together with refractive laser triangulation. Quéau et al. propose complete models for orthographic cameras and directional lights in [START_REF] Quéau | On Photometric Stereo in the Presence of a Refractive Interface[END_REF].

Inverse rendering and novel views. Inverse/differentiable rendering techniques have started to emerge recently, such as the differentiable rasterizers [START_REF] Liu | Soft rasterizer: A differentiable renderer for image-based 3d reasoning[END_REF][START_REF] Munkberg | Extracting triangular 3d models, materials, and lighting from images[END_REF], but rasterization is not well adapted to refraction. Ray-tracing/particle tracking inverse renderers are more promising [START_REF] Nimier-David | Mitsuba 2: A retargetable forward and inverse renderer[END_REF][START_REF] Jakob | JIT: a just-in-time compiler for differentiable rendering[END_REF], or [START_REF] Yan | Efficient estimation of boundary integrals for path-space differentiable rendering[END_REF] which presents a reconstruction of an object in a cuboid refractive medium. They do not seem very robust (yet) and are very computationally heavy. Adaptations of the neural radiance field framework [START_REF] Mildenhall | NeRF: Representing Scenes as Neural Radiance Fields for View Synthesis[END_REF] to refraction is proposed in [START_REF] Bemana | Eikonal fields for refractive novel-view synthesis[END_REF][START_REF] Fujitomi | Light Bending Neural Radiance Fields for Transparent Medium[END_REF], but as NeRF does, they only deal with novel view synthesis. Building on the 3D reconstruction NeuS framework by [START_REF] Wang | Neus: Learning neural implicit surfaces by volume rendering for multi-view reconstruction[END_REF], a recent article proposes an extension for an object encased in a cuboid refractive medium.

In this paper, we rather assume that the camera poses have been pre-calibrated, e.g. by structure-from-motion, and we focus on 3D reconstruction by multi-view stereo (MVS). Extensions of the PMVS algorithm (patch-based multi-view stereo) by Furukawa and Ponce [START_REF] Furukawa | Accurate, Dense, and Robust Multiview Stereopsis[END_REF] have been proposed by Kang et al. [START_REF] Kang | Two-View Underwater Structure and Motion for Cameras under Flat Refractive Interfaces[END_REF] and by Chang and Chen [START_REF] Chang | Multi-View 3D Reconstruction for Scenes under the Refractive Plane with Known Vertical Direction[END_REF], but the extension process is not really detailed. Agrawal et al. have also studied MVS in a refractive medium in [START_REF] Agrawal | A theory of multi-layer flat refractive geometry[END_REF], but their work remains limited to a (multi-layer) planar interface. To the best of our knowledge, the only work considering a non-planar interface is the recent work of Yoon et al. [START_REF] Yoon | Depth estimation from stereo cameras through a curved transparent medium[END_REF]. However, their objective is different, because the curved refractive medium is placed between the observer and the object, while in our case the object lies inside the refractive medium.

Overall, there is a lack of multi-view stereo technique adapted to the case of an object immersed in a refractive medium with a non-planar shape. The goal of our work is to fill this gap.

From Multi-view Stereo to Refracted Multi-view Stereo

Multi-view Stereo

Multi-view stereo (MVS) is a dense 3D reconstruction technique whose principle is to maximize the photometric coherence between different images in a 3D scene. An overview of the different approaches to MVS is provided in [START_REF] Furukawa | Multi-View Stereo: A Tutorial[END_REF]. Knowing t + 1 images of the scene and their associated t + 1 camera poses, we choose the image of the first pose as the reference image. For a 3D point visible in all the images, denoted P, we denote by p = π(P) its projection in the reference image and p j = π j (P), j ∈ {1, . . . , t}, its projections in the t other images, called control images. The Lambertian assumption is written:

I j • π j • π -1 z (p) p j = I(p), j ∈ {1, . . . , t} (1)
where I j and I denote the gray level functions of the j-th control image and the reference image, and where the index z of the function π -1 z is required, because the point P = π -1 z (p) is defined only if its depth z is known.

The MVS technique consists in searching for the point P = π -1 z (p), conjugate of p, satisfying the system of Equations (1), for example by solving the least squares problem:

min z ∈ R t j=1 I j • π j • π -1 z (p) -I(p) 2 (2)
In practice, the comparison between the gray levels I j and I is performed between the neighborhoods of p j and of p, and the use of a robust estimator is recommended (see the overview presented in [59, chapter 2]).

When the medium is homogeneous, the π -1 z transformation from the reference view to the 3D scene consists in inverting the central projection. Denoting by K the calibration matrix of the camera, this transformation is written:

π -1 z (p) = z K -1           p 1           (3)
The reprojection in the j-th control image is also obtained by central projection, taking into account the change in camera pose.

If the rigid transformation between the reference pose and the j-th pose is characterized by the rotation matrix R j and the translation vector t j , assumed to be known, and using a function

ρ such that ρ [a, b, c] = [a/c, b/c] ,
this second transformation is written:

π j (P) = ρ K (R j P + t j ) (4)
Carrying over ( 3) and ( 4) into ( 2), the problem of 3D reconstruction by MVS is rewritten:

min z ∈ R t j=1           I j • ρ           K (z R j K -1           p 1           + t j )           -I(p)           2 (5)
The argument to be minimized in ( 5) may be nonlinear, nondifferentiable and/or nonconvex, which makes optimization potentially difficult. For this reason, solving ( 5) is usually done by an exhaustive search (brute-force) in a predefined list of values of depth z (see Figure 2-a). This strategy, apparently simplistic, turns out to be very efficient for the 3D reconstruction of scenes with a sufficiently textured surface [START_REF] Goesele | Multi-View Stereo Revisited[END_REF]. As suggested in Figure 2-b, the scenario is similar, but more complex, when the 3D scene is immersed in a refractive medium.

p z p 1 p 2 p 3 P C C 1 C 2 C 3 p C C 1 C 2 C 3 p 1 p 2 p 3 P P 1 P 2 P 3 P (a) (b)

Refractive Multi-view Stereo

The image of an object immersed in a medium whose index of refraction (IoR) is strictly greater than 1, called refractive medium, is distorted. Among other consequences, the epipolar geometry of a pair of images is lost. Instead of corresponding to an epipolar line in the second image, a point in the first image corresponds to a curve whose shape depends on the IoR of the refractive medium and on the shape of the interface between the refractive medium and the surrounding air, called interface. Chari and Sturm have shown in [START_REF] Chari | Multi-View Geometry of the Refractive Plane[END_REF] how to generalize the matrix formalism of the epipolar geometry, using a fundamental matrix of size 12 × 12.

The epipolar geometry is mainly dedicated to the estimation of the camera pose by structure-from-motion. Since this problem has already been extended to take into account refraction [START_REF] Jordt | Refractive 3D reconstruction on underwater images[END_REF][START_REF] Chadebecq | Refractive Two-View Reconstruction for Underwater 3D Vision[END_REF], we focus in this paper on extending the MVS technique to the case where the 3D scene is immersed in a refractive medium. This new problem, called RMVS (for Refractive Multi-View Stereo), consists in solving, at each point p of the reference image, the same problem (2) as before, with several adaptations:

• Back-projection of an image point p -In the presence of refraction, the back-projection of p requires tracing the ray coming from C and passing through p, which is now a broken line (see Figure 2-b). The expression of the back-projection is thus slightly more complicated than (3):

π -1 z (p) = P + z v (6)
where the point of incidence P lies at the intersection of the deprojection line of p with the diopter, and the unit director vector v of the refracted ray is obtained by the Snell-Descartes laws of refraction (see Section 3.3). In addition, z ≥ 0 denotes the distance between P and P along the refracted ray (see Figure 2-b). The determination of the point of incidence P can be more or less tedious, depending on the shape of the interface, while the calculation of the unitary directing vector v does not present any particular difficulty, as long as the normal to the interface is known with sufficient accuracy.

• Reprojection of a 3D point P -On the other hand, computing the reprojection p j = π j (P) in the presence of refraction is considerably more complicated than (4). Instead of a simple central projection into the j-th control image, the π j function requires solving a shortest optical path problem (see Section 3.3).

• Image multiplication -As the examples in Figure 6 show, due to refraction, a single 3D point P can be projected to multiple points in the image. However, this does not really complicate the problem, since each of these projections can be used equivalently to solve [START_REF] Brument | A Shape-from-silhouette Method for 3D-reconstruction of a Convex Polyhedron[END_REF].

Compared to the MVS technique, the main difficulty of RMVS is the reprojection p j = π j (P), j ∈ {1, . . . , t}, of a 3D point P into the different control images. Let us first consider the case of a plane diopter, before tackling the case of a diopter of any shape.

Planar Interface

According to the first Snell-Descartes law on refraction, the refracted ray is in the plane of incidence, defined by the incident ray and the normal to the dioptre in P, which allows us to represent the phenomenon in a plane (see Figure 3-a).

By noting i 1 the angle between the normal to the dioptre and the incident ray, in a medium of IoR n 1 , and i 2 the angle between this same normal and the refracted ray, in a medium of IoR n 2 , the second Snell-Descartes law on refraction is written: In the case of a planar interface, squaring the two members of ( 7), and using the notations in Figure 3-a, we have:

n 1 sin i 1 = n 2 sin i 2 (7) P u 2 P C j u 1 i 1 i 2 n 2 n 1 O v 1 v 2 u v (a) (b)
n 2 1 (u 1 -u) 2 (u 1 -u) 2 + v 2 1 = n 2 2 (u 2 -u) 2 (u 2 -u) 2 + v 2 2 (8)
To determine the point of incidence P, we have to solve a quartic equation in ū:

a 4 u 4 + a 3 u 3 + a 2 u 2 + a 1 u + a 0 = 0 (9)
whose coefficients are easy to express in terms of u 1 , v 1 , u 2 , v 2 and α = n 2 /n 1 [START_REF] Agrawal | A theory of multi-layer flat refractive geometry[END_REF].

The image of a point P is unique if the interface is planar, which implies that Equation ( 9) has a unique real solution, which can be obtained e.g., by the Newton-Raphson method. The computation of p j = π j (P) is thus performed in two steps. First, determine the point of incidence P by solving [START_REF] Zhang | On-site calibration of underwater stereo vision based on light field[END_REF]. Then, project this point into the j-th control image, according to (4).

Interface of Any Shape

The Huygens-Fresnel principle allows to predict the wave surfaces of any luminous flux, which are orthogonal to the light rays.

A discrete version of the calculation of wave surfaces is provided by Dijkstra's algorithm [START_REF] Dijkstra | A note on two problems in connexion with graphs[END_REF], which allows to find the shortest path between two vertices of a graph. To compute the path of the light rays, it is therefore possible to divide the scene into voxels, which constitute the vertices of an undirected graph. However, the problem is greatly simplified if we assume that the refractive medium is homogeneous.

Indeed, if the IoR is the same at any point of the refractive medium, the light propagates in straight line, as in the air. The path of a light ray between a 3D point P and the center of projection C j , j ∈ {1, . . . , n}, is thus a broken line with only one break at the interface (see Figure 2-b). As we have already said at the end of the previous section, finding the shortest optical path between P and C j comes down to finding the point of incidence P. If the interface is planar, we know from Section 3.3 that this is equivalent to solving a quartic equation, but if the shape of the dioptre is more complex, the analytical solving is generally impossible.

It is legitimate to ask whether the π j transformation preserves the alignment of the points, i.e. whether the π j image of a refracted light ray, which is a straight line, is still a straight line in the j-th control image, in which case the conjugate point of a p point of the reference image could be searched for along this line. The example in Figure 3-b shows that this is not the case.

To find the point of incidence P in the general case, a very simple solution consists in discretizing the interface, and to minimize the optical path of the ray (C j , P, P) passing through one of the P points of the discretized interface, which amounts to solving a potentially complicated optimization problem by an exhaustive search in a list of "eligible" points P:

P = argmin P n 1 d(C j , P) + n 2 d( P, P) (10)
where d(•, •) denotes the Euclidean distance in R 3 .

In practice, the interface is discretized in the form of a 3D mesh with triangular faces, the points P eligible in the search for the point of incidence P being the barycenters of the triangles visible from the center of projection C j (see Figure 4). The triangle colored in blue on Figure 5 corresponds to the solution of Problem [START_REF] Lee | A biprism-stereo camera system[END_REF]. To refine the result, the point P is then searched in the plane defined by this triangle, according to the method described in Section 3.3. This solution is retained if it is inside the triangle.

It is rejected if it is not, and a search, identical in all points, is performed on the set of adjacent triangles, which are indicated in purple on Figure 5.

P C j P Figure 4.
The point of incidence P between a 3D point and the center of projection C j of the j-th control camera is determined by testing the set of barycenters P of the triangles of the 3D mesh of the interface that are seen by this camera.

C j P P Figure 5. Once the triangle corresponding to the solution of Problem [START_REF] Lee | A biprism-stereo camera system[END_REF] has been identified (triangle indicated in blue), the search for the point of incidence P is refined using the method described in Section 3.3. In the case where the solution of this second problem is outside the triangle, a search is performed on the set of adjacent triangles (triangles indicated in purple).

Validation on Synthetic Images

Cubic Interface

We start by validating the proposed method on a scene consisting of an insect, called graphosoma, with a size of about 30 mm, immersed in a refractive cube whose IoR is equal to that of the epoxy resin (n 2 = 1.56). The focal length of the camera is equal to 50 mm, and the average distance from the scene to the camera is about 180 mm. Figure 6 shows two synthetic images (among 18) of this scene, computed with the ray tracing functions of the software Blender. As a comparison, Figure 9 shows the result obtained by MVS from the same data, i.e. without taking refraction into account: because of the demultiplication of the image of the graphosoma, it is badly reconstructed, and three copies even appear. On the other hand, this improvement is obtained at the cost of a much higher computation time: we go from 16 minutes for the result of Figure 9, to 24 hours for the one of Figure 7 (CPU Intel Xeon Silver 4110 2.10 GHz using all 32 threads on parallel computing), for a number of 3D points of the order of 500.000 in both cases. Note however that the computation time could be greatly reduced if, among the barycenters P of the mesh triangles (see Figure 4), only those that project inside the insect silhouette, in the considered control image, were tested.

This first example allows us to detail the solving method. The 3D reconstruction of Figure 7 results from the merging of eight colored 3D point clouds, each of these clouds being obtained as follows: • One of the images is chosen as reference image. Five other images are chosen as control images: for four of them, the "main" image of the insect is seen through the same face of the cube as in the reference image; for the fifth control image, it is an adjacent face.

• For each pixel p of the reference image, for each point P located on the refracted ray resulting from the deprojection of p, for each control image, and for each face of the cube visible in the control image (there are at most three faces visible per control image), a quartic equation of type ( 9) is solved by the Newton-Raphson method.

• Once the solution has been found, if its projection p j in the j-th control image is located inside the silhouette of the insect in the face concerned, then the similarity between the neighborhoods of p and p j is calculated, using a sufficiently robust estimator, in this case the sum of absolute deviations (SAD).

• In the case where the SAD can be calculated for several faces of the j-th control image (at most three), only the smallest value is retained.

• Otherwise, that is to say if no SAD can be calculated, another point P is tested, among a predefined list of 3D points.

• The retained point P is the one which minimizes the SAD. Its color is that of the pixel p of the reference image.

• If no SAD can be calculated, no 3D point is associated with pixel p.

Given that all the images of the graphosoma are synthetic, it is possible to measure the deviations from the ground truth of each of the eight 3D point clouds whose fusion produces the result of Figure 7. The square root of the mean and the square root of the median of the squared deviations are reported in Table 1. These values are very low, compared to the scale of the reference 3D model and the distance from the 3D model to the camera, which allows us to quantitatively validate our RMVS solving method. Figure 10-a shows that the image of a point P located inside the refractive medium can be multiplied. These multiple images correspond to as many local minima of the optical path between a 3D point P and the center of projection (Fermat's principle). It is therefore possible to match these different images of the insect, which amounts to applying our RMVS solving method using only one view, as was done in [START_REF] Lee | A biprism-stereo camera system[END_REF][START_REF] Yamashita | Monocular Underwater Stereo -3D Measurement Using Difference of Appearance Depending on Optical Paths[END_REF]. The result (see

Face

Spherical Interface

The second test is carried out on a graphosoma immersed in a sphere of epoxy resin. Figure 11 shows two synthetic images of this new scene, and puts them next to images of the graphosoma seen from the same two angles, but outside the refractive medium.

Besides the magnifying effect due to the convex shape of the block of epoxy resin, which acts as a convex lens, the main difference between these two pairs of images lies in the legs and insect antennae, which seem deformed by refraction. It is noteworthy that the images in Figures 11-a For comparison, Figure 13 shows that, without taking refraction into account, MVS fails to reconstruct the 3D shape, even roughly. Ghosted legs and antennae appear, which reflects the lack of consistency between the eight 3D point clouds. The interface discretization scale is guided by the need to find a compromise between precision and computing time. Table 2 shows the influence of a decrease in the number of triangular faces of the 3D mesh of the sphere, which means a less precise knowledge of the interface, through the two estimators already used in Section 4. Table 2. Influence of a decrease in the number of triangular faces of the 3D mesh of the sphere on the 3D reconstruction of the graphosoma by our RMVS solving method (the numbers of faces are indicated in parentheses).

Figure 14 shows the 3D reconstructions corresponding to Table 2. As expected, the 3D points which "disappear" first, i.e., those conjugated with pixels p for which no SAD similarity value can be calculated, are located on the thinnest parts of the 3D model, namely the legs and the antennae. 2). Figure 12 shows the result when all the 327,000 faces are used.

Other Interfaces

Figure 15 shows two synthetic images of the graphosoma immersed in a regular dodecahedron made of the same epoxy resin.

The 3D reconstruction obtained by our RMVS solving method is illustrated in Figure 16. The deviations from the ground truth are hardly higher than those of Table 1, even though the images used are much more difficult to interpret, in terms of 3D, than those of Figure 7: if the RMSE, computed on the whole point cloud, goes from 0.57 mm to 1.10 mm, the RMedSE increases much less, from 0.30 mm to 0.35 mm. This increased value of the RMSE probably comes from the significant deformation of the images, so that the similarity measurement by the SAD estimator is somehow biased. Figure 17 shows that images can be distorted even more with a block of any convex shape. The 3D reconstruction obtained by our RMVS solving method (see Figure 18), which remains faithful to the original, is still less precise than that of Figure 7, because for some grazing rays (see Figure 3-a), the angle i 2 of the Snell-Descartes law ( 7) is very close to π/2. Hence, as the derivative of function arcsin tends to infinity at 1, this leads to inaccuracies in the calculation of the angle i 1 of Equation [START_REF] Agrafiotis | DepthLearn: learning to correct the refraction on point clouds derived from aerial imagery for accurate dense shallow water bathymetry based on SVMs-fusion with LiDAR point clouds[END_REF]. By comparison, In this section, we validated our RMVS solving method only on synthetic images, but purposely. Indeed, to be able to process real images, several additional data are necessary, in addition to the images themselves and the intrinsic parameters of the camera:

the shape of the interface (with more or less precision, see Table 2), the poses of the camera and the IoR of the transparent medium.

Implementation on Real Images

The main difficulty in implementing our 3D reconstruction method on real images is estimating the camera poses. We could use the refractive structure-from-motion method proposed in [START_REF] Chadebecq | Refractive Two-View Reconstruction for Underwater 3D Vision[END_REF], but this would require knowing the IoR of the medium, which is one of the unknowns. As we also need to recover the interface 3D shape, we propose in this section a method of simultaneously estimating the poses of the camera and the interface 3D shape. The method combines shape-from-silhouette [START_REF] Hernández | Silhouette and Stereo Fusion for 3D Object Modeling[END_REF] with point correspondences across the views and does not involve the IoR. The latter can therefore be estimated a posteriori, as we will see in Section 5.2.

Estimating the Camera Poses and the Interface 3D Shape

In this subsection, we describe a method for collecting the camera poses and the 3D shape of the interface, expressed in a common 3D frame, consisting of fixing the camera in front of the object to be analyzed, positioned on a turning table. The rotation of the table acts as if the camera was rotating around the object. The simultaneous estimation of the poses of the camera and the interface 3D shape is performed by shape-from-silhouettes.

This technique has the advantage of not being dependent on the IoR of the refractive medium. Nevertheless, it is hopeless to solve the problem for an object of any shape. Indeed, the principle of shape-from-silhouettes is to compute a volume enclosing the object by intersecting the deprojections of its silhouettes. However, to hope for an accuracy that meets our needs, it would be necessary to have an infinite number of poses, unless we restrict ourselves to the particular case of a polyhedral interface, whose shape can be characterized by a small number of vertices. In the case of insects trapped in amber, this is quite realistic since it is possible, as already said in Section 1, to cut the amber in order to give it a polyhedron shape i.e., a a surface composed of several planar faces. Now, we consider a scene consisting of a convex polyhedron with N vertices, placed on a turntable. The vertices coordinates X n , n ∈ {1, . . . , N}, are expressed in a 3D frame R ref attached to the turntable. The origin of R ref is located at the intersection between the supporting plane of the table and the rotation axis, and its two first axes define the supporting plane of the table, such that the third axis is a normal vector to the table directed upwards. A series of J views of this scene is acquired by a static camera of known intrinsics. In the j-th image, j ∈ {1, . . . , J}, the image x j n of vertex X n satisfies the equation:

x j n ∼ P j X n [START_REF] Yamashita | Monocular Underwater Stereo -3D Measurement Using Difference of Appearance Depending on Optical Paths[END_REF] where the perspective projection matrix P j related to the j-th view is written:

P j = K [ R ref R θ j | t ref ] (12)
In [START_REF] Chen | Depth from refraction using a transparent medium with unknown pose and refractive index[END_REF], K is the calibration matrix, (R ref , t ref ) define the change of 3D coordinates from R ref to the camera frame, and matrix R θ j denotes the rotation of the table by an angle θ j around its rotation axis:

R θ j =                  cos θ j -sin θ j 0 sin θ j cos θ j 0 0 0 1                  (13)
In homogeneous Cartesian coordinates x = [x, y, 1] , an ellipse can be represented as x A x = 0, where A is a symmetric 3 × 3 matrix under suitable conditions on its coefficients. The interior (resp. exterior) of the ellipse is given by the points x for which

x A x < 0 (resp. x A x > 0). A first step consists in creating polygonal silhouettes of each view. This is obtained by simple operations: background subtraction from a reference image, thresholding, morphological processing, extraction, and simplification of the convex hull. Then the collection V of all the polygon vertices for all the silhouettes is created. In a second step, V must be robustly partitioned in subsets of points located on common ellipses, as these ellipses should be the images of parallel circular trajectories in 3D space. The partition size, which should match the number of vertices of the polyhedron, is unknown. A solution to this partitioning problem is described in [START_REF] Brument | A Shape-from-silhouette Method for 3D-reconstruction of a Convex Polyhedron[END_REF]. The authors propose to exploit the parallelism of the vertices trajectories. The images of the turntable circular points (ICP) [START_REF] Hartley | Multiple View Geometry in Computer Vision[END_REF], noted h 1 ± ih 2 , are estimated along with the correspondences and are assumed to be known is the sequel. See [START_REF] Brument | A Shape-from-silhouette Method for 3D-reconstruction of a Convex Polyhedron[END_REF] for more details.

Assuming that the ICP h 1 ± ih 2 , the correspondences {x j n } and the calibration matrix K are known, the problem consists in estimating the positions of the vertices and the poses of the polyhedron in the camera frame. More precisely, it is to compute the rotation matrix R ref , the translation vector t ref , the 3D coordinates of the vertices {X n } n∈{1,...,N} and the angles {θ j } j∈{1,...,J} . Computing R ref and t ref is performed using the method proposed in [START_REF] Sturm | Algorithms for Plane-Based Pose Estimation[END_REF].

The angle of rotation θ j is the angle of rotation of the turntable from a reference position, for example at the time of the first acquisition, and its value in view number j. Angle θ j is computed as:

θ j = j-1 k=1 θ k,k+1 (14)
where θ k,k+1 is the angle of rotation between two consecutive acquisitions k and k + 1. The value of θ k,k+1 is obtained by computing the median of the cosines of the estimated angles from the visible vertices:

θ k,k+1 = acos median l∈D k d k l d k+1 l (15)
where D k ⊂ {1, . . . , N} denotes the set of vertex indices that have been detected in the k and k + 1 images, and where d k l and d k+1 l denote the unit vectors pointing towards the images by H -1 of the corresponding points x k l and x k+1 l in the consecutive views k and k + 1, that is:

d k l = f (H -1 x k l ) -f (H -1 c l ) f (H -1 x k l ) -f (H -1 c l ) (16)
and likewise for d k+1 l . In ( 16), function f is such that f ([u, v, w] ) = [u/w, v/w] , and c l is the homogeneous coordinate vector of the image of the centre of the trajectory, supposedly circular, of the vertex number l, obtained from the pole-polar relationship

c l = A -1 l l ∞
, where l ∞ denotes the vanishing line vector of the table plane given by the cross-product l ∞ = h 1 × h 2 , and A l the matrix of the ellipse image of the trajectory of the vertex number l. The table is assumed to rotate in the trigonometric direction during the acquisition. The values of θ k,k+1 are supposed to be between 0 and 180 degrees.

At this stage, all the camera poses are known. The 3D coordinates of the vertices {X n } n∈{1,...,N} are obtained by triangulating the correspondences {x j n }. Camera poses and vertices coordinates are then refined through a bundle adjustment minimizing the Euclidean distances between the correspondences and their reprojections.

Validation on Real Images

The estimation of the index of refraction (IoR) can be performed using a dedicated instrument called refractometer, but we propose in this subsection an estimation method that takes advantage of the joint estimation of the camera poses and of the interface 3D shape, described in the previous subsection. Indeed, our RMVS solving method, which is described in Section 3 and validated on synthetic images in Section 4, can be tested for different values of the IoR, for example in a range around its "plausible" value. It remains to find a sufficiently discriminating criterion for this estimation. Figure 22 shows that the number of effectively reconstructed 3D points constitutes such a criterion.

We are now able to test the complete RMVS resolution pipeline on real data, on the sole condition that the interface is polyhedral.

Figures 23 and24 show two objects on which tests were carried out. These are two parallelepipeds made of resin containing, respectively, a beetle and a grasshopper. Figure 25 shows how the previous method, applied to the images of the beetle, gives an estimate of the resin IoR equal to n 2 = 1.50. The 3D reconstructions of these two insects are presented in Figures 26 and27. Clearly, the beetle is better reconstructed than the grasshopper. This is mainly due to the characteristics of the resin block containing the beetle, which do not perfectly satisfy the assumptions on which our RMVS solving method is based. On the one hand, at least one of the faces is not as planar as the others. On the other hand, the resin consists of layers visible to the naked eye, which suggests that the light rays inside the refractive medium are not perfectly straight. But another important difference between the results of Figures 26 and 27 concerns the insects themselves. Indeed, some parts of the grasshopper's body appear relatively translucent, which calls into question one of the basic assumptions of the MVS technique and its variants, namely that the surface to be reconstructed is assumed to be opaque and Lambertian.

A final test carried out on the images of the beetle aims to verify, at least qualitatively, to what extent the use of an erroneous IoR value influences the result. Figure 28 shows two 3D reconstructions of the beetle, obtained with two different values of the IoR:

the result on the right, obtained with a slightly overvalued IoR (n 2 = 1.56), is clearly less precise than the one on the left, where the value of the IoR (n 2 = 1.50) was estimated using the method described above (see Figure 22).

Conclusion and Perspectives

In this paper, we have shown how to adapt the multi-view stereo technique to the case where the object of interest is immersed in a refractive medium. Knowing that refraction distorts the images, it is necessary to take this phenomenon into account to model the path of the light rays. More specifically, we have proposed a fully discrete RMVS solving method. Our first results on real data are encouraging, although there are still many obstacles to overcome before this study leads to an operational method that can be used by entomologists.

As has already been said, the relatively limited quality of the result of Figure 27 comes from both the block of resin and the object it contains, neither of which fully satisfies the assumptions on which our RMVS solving method is based. A first perspective is therefore to make the entire pipeline more robust to the most predictable defaults. In particular, the problem becomes much more complex if the index of refraction of the medium is not homogeneous. Moreover, even under the Lambertian assumption, if the refractive medium is colored, the appearance of a 3D point can vary from one image to another, depending on the distance traveled by the light in the medium. Finally, focus blur, which is a phenomenon difficult to avoid at small scales, was ignored in this study.

Taking all these effects into account will certainly improve the accuracy of the results.

The efficiency of our RMVS solving method could also be improved. To give an order of magnitude of the complexity of this fully discrete solving strategy, each of the eight 3D point clouds which give, by fusion, the result of Figure 7, is obtained with t = 5 control images. The reference image has approximately 1.5 × 10 5 pixels. The block of epoxy resin being discretized in 200 slices, the complete resolution of the problem (2) requires to solve approximately 10 8 equations of degree 4, which requires around 24 hours of CPU time.

Figure 1 .

 1 Figure 1. Left: prehistoric beetle trapped in amber (seen under a microscope). Right: reptiles specimens in jars (image from the Natural History Museum in Copenhagen).

Figure 2 .

 2 Figure 2. (a) MVS in a homogeneous medium: the different proposals for the point P, which are materialized by red dots, are reprojected in the control images. (b) MVS with refraction: the reprojection of P in the control images is more difficult to compute, due to refraction.

Figure 3 .

 3 Figure 3. (a) Second Snell-Descartes law on refraction. (b) Synthetic images of two straight lines (in red and in blue) crossing a refractive sphere.

Figure 6 .

 6 Figure 6. Two synthetic images (among 18) of a graphosoma immersed in a cube of epoxy resin. In both cases, the insect is visible through three faces of the cube (only partially, regarding the top face). The other fragments of images of the insect visible near the edge of the cube silhouette, which are due to reflection phenomena, are not used by our solving method. Source of the 3D model: Digital Archive of Natural History [62].

Figure 7

 7 Figure7shows three views of the colored 3D point cloud reconstructed by our RMVS solving method. Figure8shows the

Figure 7 .

 7 Figure 7. 3D reconstruction of a graphosoma immersed in a cube of epoxy resin, seen from three angles, obtained by our RMVS solving method from 18 synthetic images such as those in Figure6.

Figure 8 .

 8 Figure 8. Result of Figure 7 after "cleaning" the 3D point cloud by the Connected-component labeling tool of the Cloud Compare software.

Figure 9 .

 9 Figure 9. 3D reconstruction obtained by MVS, without taking refraction into account: because of the demultiplication of the image, the graphosoma is badly reconstructed, and it even appears copied three times.

Figure 10 .

 10 Figure 10. (a) Example where the image of the graphosoma is detripled. A 3D point cloud can be calculated from this single image, by choosing the "main" image (right face) as reference image. (b-c) Two views of the 3D point cloud reconstructed by our RMVS solvind method, from this single image.

Figure 11 .

 11 Figure 11. (a-b) Two synthetic images of the graphosoma immersed in a sphere of epoxy resin.(c-d) Synthetic images of the graphosoma seen from the same two angles, but outside the refractive medium. In addition to the magnifying effect due to the convex shape of the block of epoxy resin, legs and antennae appear deformed.

Figure 12 shows

 12 Figure12shows from three angles the colored 3D point cloud obtained by our RMVS solving method, from 18 images such

Figure 12 .

 12 Figure 12. 3D reconstruction of the graphosoma immersed in a sphere of epoxy resin, seen from three angles, obtained by our RMVS solving method, from 18 images such as those in Figures 11-a and 11-b, after cleaning by the Connected-component labeling tool of the Cloud Compare software.

Figure 13 .

 13 Figure 13. 3D reconstruction by MVS, from 18 images such as those in Figures 11-a and 11-b: due to the lack of consistency between the eight 3D point clouds, ghost legs and antennae appear.

Figure 14 .

 14 Figure 14. Six 3D reconstructions of the graphosoma immersed in a sphere of epoxy resin, when the percentage of triangular faces of the sphere used by our RMVS solving method decreases (see Table2). Figure12shows the result when all the 327,000 faces are used.

Figure 15 .

 15 Figure 15. Two images of the graphosoma immersed in a regular dodecahedron of epoxy resin.

Figure 16 .

 16 Figure 16. 3D reconstruction of the graphosoma immersed in a regular dodecahedron of epoxy resin, seen from three angles, obtained by our RMVS solving method from 18 images such as those in Figure 15, after cleaning with the Connected-component labeling tool of the Cloud Compare software.

Figure 17 .

 17 Figure 17. Two images of the graphosoma immersed in a convex block of epoxy resin.

Figure 18 .

 18 Figure 18. 3D reconstruction of the graphosoma immersed in a convex block of epoxy resin, seen from three angles, obtained by our RMVS solving method from 18 images such as those in Figure17.

Figure 19 .

 19 Figure 19. 3D reconstruction by MVS, from 18 images such as those in Figure 17: the result looks like a random 3D point cloud.

Figure 20 .

 20 Figure 20. (a) Overview of the acquisition setup. The polyhedron is placed on a turntable. The views are acquired by a static camera. The origin t ref is the intersection of the table with its axis of rotation. The rotation matrix R ref , with columns r 1 , r 2 and r 3 , describes the pose of the turntable with respect to the camera frame. (b) The image of a vertex, located at a distance ρ n of the rotation axis, describes an ellipse A n in the pixel plane. Its imaged center c n satisfies the pole-polar relationship c n = A -1 n l ∞ .

Figure 21 .

 21 Figure 21. (a) One out of 40 images of a parcel placed on a turntable. The silhouettes of the parcel, considered as a convex polyhedron, are extracted in all views. (b) All the silhouettes vertices, displayed in red.

Figure 22 .

 22 Figure 22. Evolution of the percentage of reconstructed 3D points in function of the index of refraction of the medium in which the graphosoma is immersed.The IoR used when simulating the images (n 2 = 1.56) corresponds precisely to the maximum of this curve, which validates the proposed criterion.

Figure 23 .

 23 Figure 23. Left: two real images of a beetle immersed in a parallelepipedic block of resin, placed on a turntable. Right: zooms on the block.

Figure 24 .

 24 Figure 24. Left: two real images of a grasshoper immersed in a parallelepipedic block of resin, placed on a turntable. Right: zooms on the block.

Figure 25 .

 25 Figure 25. Evolution of the percentage of reconstructed 3D points, in function on the IoR of the medium in which the graphosoma is immersed. The maximum of this curve gives us an estimate of the IoR equal to n 2 = 1.56.

Figure 26 .

 26 Figure 26. 3D reconstruction of the beetle from 24 images such as those in Figure23, using our RMVS solving method.

Figure 27 .

 27 Figure 27. 3D reconstruction of the grasshoper from 24 images such as those in Figure24, using our RMVS solving method.

Figure 28 .

 28 Figure 28. Comparison of our RMVS solving method, tested on the same 24 real images of the beetle (see Figure23), using either the value of the IoR estimated by the method illustrated in Figure22(n 2 = 1.50), or a slightly overvalued IoR (n 2 = 1.56). The first result is obviously more accurate.

Table 1 .

 1 Second line: root mean square error (RMSE, in mm) of the eight 3D point clouds whose fusion provides the result of Figure7. Third line: root median square error (RMedSE, in mm). The last column shows the average of these estimates across the eight 3D point clouds.

		Front	Front-right	Right	Back-right	Back	Back-left	Left	Front-left	All
	RMSE (mm)	0.25	0.48	0.85	0.98	0.40	0.73	0.60	0.65	0.57
	RMedSE (mm)	0.15	0.30	0.40	0.40	0.25	0.33	0.33	0.40	0.30

  1 (RMSE and RMedSE), as well as the percentage of 3D points actually reconstructed, and the CPU time.

	Percentage of faces	100 (327k)	50 (164k)	25 (82k)	10 (33k)	5 (16k)	1 (1.6k)	0.5 (820)
	RMSE (mm)	1.77	1.82	1.90	2.03	2.25	3.72	4.32
	RMedSE (mm)	0.30	0.33	0.35	0.45	0.57	2.10	3.05
	Percentage of 3D points	98.0	97.4	96.2	93.9	91.6	77.8	69.4
	CPU time (minutes)	498	272	147	49	37	28	23