Riadh Ben Salah
email: riadhbs@yahoo.com

Olivier Alata
email: olivier.alata@univ-st-etienne.fr

Benoit Tremblais
email: benoit.tremblais@univ-poitiers.fr

Lionel Thomas
email: lionel.thomas@univ-poitiers.fr

Laurent David
email: laurent.david@univ-poitiers.fr

Tomographic Reconstruction of 3D Objects using Marked Point Process Framework

Keywords: 3D Object Reconstruction, Stochastic Models, Marked Point Processes, Metropolis-Hastings-Green (MHG) algorithm, Simulated annealing, Tomographic PIV

For reconstructing sparse volumes of 3D objects from projection images taken from different viewing directions, several volumetric reconstruction techniques are available. Most popular volume reconstruction methods are algebraic algorithms (e.g. the multiplicative algebraic reconstruction technique, MART). These methods which belong to voxel-oriented class, allows volume to be reconstructed by computing each voxel intensity. A new class of tomographic reconstruction methods, called "object-oriented" approach, has recently emerged and was used in the Tomographic Particle Image Velocimetry technique (Tomo-PIV). In this paper, we propose an object-oriented approach, called IOD-OVRMPP (Iterative Object Detection -Object Volume Reconstruction based on Marked Point Process), to reconstruct the volume of 3D objects from projection images of 2D objects. Our approach allows the

problem to be solved in a parsimonious way by minimizing an energy function based on a least squares criterion. Each object belonging to 2D or 3D space is identified by its continuous position and a set of features (marks). In order to optimize the population of objects, we use a simulated annealing algorithm which provides a "Maximum A Posteriori" (MAP) estimation. To test our approach, we apply it to the field of Tomo-PIV where the volume reconstruction process is one of the most important steps in the analysis of volumetric flow. Finally, using synthetic data, we show that the proposed approach is able to reconstruct densely seeded flows.

Introduction and Motivation

In the 80s, tomographic reconstruction methods were generally classified into two categories: analytical methods [START_REF] Censor | Finite series-expansion reconstruction methods[END_REF] and algebraic methods [START_REF] Lewitt | Reconstruction algorithms: transform methods[END_REF]. Analytical methods are fast and simple to implement. They give high quality reconstruction results only if the trajectory acquisition is full with a good sampling of projection images. However, with a limited number of views or noisy images, these methods perform poorly. In this case, it is better to use algebraic methods which adopt a discrete approach to the problem. Algebraic methods assume that the physical volume to reconstruct is sampled on a regular 3d grid. Thus, the volume is formed of voxels with unknown values. The assumption is that each detected value in a pixel from projection images is a linear combination of voxel values. This physically means that each voxel contributes (possibly with a zero coefficient) to each line of sight.

The estimation of the volume intersection can be readily performed assuming voxels as cubes or spheres. In this formulation, the problem involves solving the system P = Rf , where P is a vector containing the set of projection images, f is a vector containing the voxel values of the volume to reconstruct and R is a matrix of coefficients of intersection that links P to f . An attemp could be made to solve the system using the generalized inverse of R: f = (R t R)

-1 R t P . However, this inverse problem admits no solution because P is noisy (measurement noise, gap between the model and the data, etc.). The resolution of this "ill-posed problem" [START_REF] Hadamard | Le probleme de Cauchy et les équations aux dérivées partielles linéaires hyperboliques[END_REF] can be achieved by using iterative tomographic reconstruction methods like algebraic approaches. Algebraic methods were introduced by Kaczmarz [START_REF] Kaczmarz | Angenäherte auflösung von systemen linearer gleichungen[END_REF]. The leader of this family of algorithms is the ART method (Algebraic Reconstruction Technique) [START_REF] Gordon | Algebraic reconstruction techniques (ART) for three-dimensional electron microscopy and X-ray photography[END_REF]. Many other algebraic methods were proposed like MART (Multiplicative ART) [START_REF] Herman | Iterative reconstruction algorithms[END_REF], AART (Adaptive ART) [START_REF] Lu | Adaptive algebraic reconstruction technique[END_REF], SIRT (Simultaneous Iterative Reconstruction Technique) [START_REF] Benson | Modified simultaneous iterative reconstruction technique for faster parallel computation[END_REF], SART (Simultaneous ART) [START_REF] Andersen | Simultaneous algebraic reconstruction technique (SART): a superior implementation of the ART algorithms[END_REF] and SMART (Simultaneous MART) [START_REF] Darroch | Generalized iterative scaling for log-linear models[END_REF][START_REF] Schmidlin | Iterative separation of sections in tomographic scintigrams[END_REF]. Later, methods based on reconstruction by blocks, such as BIMART (Block Iterative MART), were introduced by Charles Byrne [START_REF] Byrne | Block-iterative algorithms[END_REF].

Algebraic methods allow the introduction of some constraints on the solution like the geometry of the acquisition system. Unfortunately, they converge, most of the time, to least squares solutions, which suppose Gaussian noise on the data. To tackle these problems, statistical methods were proposed. They consider P , and sometimes f , as random variables, to which a probability distribution reflecting uncertainty in the data (noise, modeling error) and f (in the case of Bayesian methods) are assigned [START_REF] Peyrin | Introduction à l'imagerie tomographique 2D et 3D reposant sur une propagation en ligne droite. cas de la tomographie par rayon X, par émission et par ultrasons[END_REF]. The best known statistical method is called MLEM (Maximum Likelihood Expectation Maximization) and was developed in 1982 by Shepp and Vardi [START_REF] Shepp | Maximum likelihood reconstruction for emission tomography[END_REF]. The MLEM method is a special case of the more general approach called Expectation Maximization (EM) [START_REF] Dempster | Maximum likelihood from incomplete data via the EM algorithms[END_REF]. An accelerated version of MLEM, called OSEM (Ordered Subsets Expectation Maximization) [START_REF] Hudson | Accelerated image reconstruction using ordered subsets of projection data[END_REF], was also proposed. Many other statistical methods were also developed like GC (Conjugate Gradient) [START_REF] Lascaux | Analyse numérique matricielle appliquée à l'art de l'ingénieur[END_REF], MENT (Maximum Entropy) [START_REF] Minerbo | MENT: A maximum entropy algorithm for reconstructing a source from projection data[END_REF] and MEM (Maximum Entropy on the Mean) [START_REF] Rietch | The maximum entropy approach to inverse problems[END_REF][START_REF] Gamboa | Méthode du maximum d'entropie sur la moyenne et applications[END_REF]. Both statistical and algebraic methods are discrete tomographic reconstruction methods.

Recently, existent tomographic reconstruction methods were applied in the fluid mechanics domain to reconstruct 3D particle volumes from a set of particle images. The images are acquired with a 3D tomographic optical system called Tomographic Particle Image Velocimetry (Tomo-PIV) which was introduced by Elsinga et al [START_REF] Elsinga | Tomographic particle image velocimetry[END_REF].

The principle of the Tomo-PIV is based on the computation of the velocity vector field of a flow from the displacement of the particles recorded on several images. Particles are seeded in the flow and are illuminated by a thick laser sheet which defines the volume of interest. The light scattered by the particles is captured by several high-resolution digital cameras from different viewing angles (usually 4 -see Fig. 1). The information on the line of sight of each pixel in a camera, through out the examined volume, is described by a polynomial approximation made from a 3D calibration procedure. Thus, the projection images set captured at a given time is supplied to a 3D tomographic reconstruction algorithm to reconstruct a volume of 3D particles. Finally, the sequences of 3D particle volumes are used to estimate the sequences of velocity fields using advanced cross-correlation techniques. The problem of tomographic reconstruction in the Tomo-PIV field is different than in other domains for the following reasons:

-The volume of interest is composed of 3D particles with a high density level (up to 800,000 particles); -The presence of ghost particles (artifacts) in reconstructed volumes due to the limited number of views (usually 4 cameras); -The experimental projection image size can reach 4000 × 4000 pixels (volume size of 4000 × 4000 × 400 voxels for example); -Noisy data due to the complex installation and vibration of the acquisition system.

In such a context, the discrete tomographic reconstruction methods proposed to solve the specific problem of TomoPIV can be divided into two classes. The first class is called fast methods and it covers algorithms designed to reduce the computational costs and accelerate the convergence to a solution by estimating the initial volume intensity distribution. The study conducted by Worth and Nickels [START_REF] Worth | Acceleration of tomo-piv by estimating the initial volume intensity distribution[END_REF] which yielded an approach named MFG (Multiplicative First Guess) can be mentioned. Atkinson and Soria [START_REF] Atkinson | An efficient simultaneous reconstruction technique for tomographic particle image velocimetry[END_REF] developed the MLOS (Multiplicative Line Of Sight) method to accelerate the MART method which is the most widely used algebraic technique in the Tomo-PIV field [START_REF] Dean | Limited-data computed tomography algorithms for the physical sciences[END_REF][START_REF] Thomas | Optimization of the volume reconstruction for classical tomo-piv algorithms (mart, bimart and smart): synthetic and experimental studies[END_REF]. Maas & al. [42] proposed the MinLOS (Minimum Line Of Sight) method which is a variant of MLOS. In the same context, other methods were developed like SAPIV (Synthetic Aperture Particle Image Velocimetry) [START_REF] Belden | Three-dimensional synthetic aperture particle image velocimetry[END_REF], MG (Multi-Grid) [START_REF] Discetti | Fast multi-resolution 3D PIV with direct correlations and sparse arrays[END_REF] and VS-TPIV (Volume Segmentation -Tomographic Particle Image Velocimetry) [START_REF] Ziskin | Volume segmentation tomographic particle image velocimetry[END_REF]. The second class is based on sparse models. In this class, the problem of tomographic reconstruction was treated in a framework exploiting the sparsity of the volume as in methods based on the L 1 norm [START_REF] Donoho | Counting the faces of randomly-projected hypercubes and orthants, with applications[END_REF] and others based on basis-pursuit algorithms [START_REF] Petra | 3D tomography from few projections in experimental fluid dynamics[END_REF][START_REF] Herzet | Bayesian pursuit algorithms[END_REF].

Analytical and discrete tomographic reconstruction methods called "voxel oriented" methods are based on the calculation of the intensity of each voxel in the volume [START_REF] Ben-Salah | Élaboration d'une méthode tomographique de reconstruction 3D en vélocimétrie par image de particules basée sur les processus ponctuels marqués[END_REF]. They take up a lot in memory storage and do not take sufficiently into account the particular shape of the objects to be reconstructed.

To tackle these problems, "object-oriented" methods, have emerged [START_REF] Wieneke | Iterative reconstruction of volumetric particle distribution[END_REF][START_REF] Champagnat | Tomographic piv: particles versus blobs[END_REF][START_REF] Ben-Salah | Élaboration d'une méthode tomographique de reconstruction 3D en vélocimétrie par image de particules basée sur les processus ponctuels marqués[END_REF]. Methods belonging to this category were developed in the fluid mechanics domain. The proposed solutions are based on the parsimony of volumes: they consider the volume as a set of 3D objects. The objects represent 3D particles which are identified by their positions and features. In 2012, Wieneke [START_REF] Wieneke | Iterative reconstruction of volumetric particle distribution[END_REF] proposed an algorithm called IPR (Iterative Particle Reconstruction) to reconstruct 3D-particle locations by comparing the recorded images with generated ones calculated from the particle distribution in the volume. In 2014, Champagnat et al. proposed the PVR (Particle Volume Reconstruction) method [START_REF] Champagnat | Tomographic piv: particles versus blobs[END_REF]. This approach is a combination between an "object-oriented" approach and algebraic methods (MART, SMART). The results obtained with numerical simulations showed that the proposed methods give equivalent reconstruction quality to the MART method for very high particle seeding density. However, proposed "object-oriented" methods are very specific to Tomo-PIV problems and cannot be exploited in other fields using sparse volumes.

In this paper, we present a general "object-oriented" approach to reconstruct volumes of 3D objects. The proposed approach can be used in several areas and meets the following criteria:

-It exploits the parsimony of the volumes using an appropriate mathematical formalism; -It allows a priori information to be easily introduced; -It works with real-space coordinates in order to locate the objects with high accuracy; -It solves an ill-posed problem even in the case of very high densities of objects; -It improves the reconstruction quality with respect to state-of-the-art ART methods and reduces memory usage as well as the artifact rate by using a formalism suited to the problem addressed.

This paper focuses on the reconstruction of volumes of 3D objects in large configuration spaces, formulated in a marked point process reconstruction framework.

Stochastic models based on marked point processes have been used in various application domains, from 3D reconstruction to texture modeling such as the reconstruction of buildings [START_REF] Ortner | Building outline extraction from digital elevation models using marked point processes[END_REF] or vascular trees [START_REF] Sun | Marked point process for vascular tree extraction on angiogram[END_REF], and the detection of tree crowns [START_REF] Perrin | Adaptive simulated annealing for energy minimization problem in a marked point process application[END_REF], populations of birds [START_REF] Descamps | Automatic flamingo detection using a multiple birth and death process[END_REF] or road networks [START_REF] Tournaire | Rectangular road marking detection with marked point processes[END_REF]. These models have proved their efficiency for object extraction in large sample spaces.

Considering previously published papers [START_REF] Ben-Salah | 3d particle volume tomographic reconstruction based on marked point process: Application to tomo-piv in fluid mechanics[END_REF][START_REF] Ben-Salah | Particle volume reconstruction based on a marked point process and application to tomo-piv[END_REF], it provides additional details and newer contributions:

-We introduce a general framework in order to make geometrical object reconstruction based on marked point processes using a minimum mean square error estimation criterion. -As one of the main objectives of Tomo-PIV is to reconstruct a velocity field using two consecutive volumes, a comparison is made on this topic between the proposed method and an optimized MinLOS-MART (Minimum Line-Of-Sight, Multiplicative Algebraic Reconstruction Technique) algorithm [START_REF] Thomas | Optimization of the volume reconstruction for classical tomo-piv algorithms (mart, bimart and smart): synthetic and experimental studies[END_REF] for synthetic data.

It should be noted that we have not provided our results on experimental data but these can be found in [START_REF] Ben-Salah | Reconstruction de volumes de particules par processus ponctuels marqués[END_REF].

This paper is organized as follows: in section 2, after a quick introduction on marked point processes, we present the general formalism of the reconstruction method, called IOD-OVRMPP (Iterative Object Detection -Object Volume Reconstruction based on Marked Point Process) and we show how to simulate it. In section 3, we present an application of our approach to reconstruct a population of 3D particles. We compare it to the reference algebraic reconstruction technique, MinLOS-MART, using synthetic data with different seeding densities. Finally in section 4, the ma-jor findings are summarized, some concluding remarks made and some perspectives given.

2 Object reconstruction based on marked point processes

Definitions

In this section, the basic ideas of point processes (PP) and marked PP (MPP) are recalled [START_REF] Van-Lieshout | Markov point processes and their applications[END_REF][START_REF] Daley | An Introduction to the Theory of Point Processes: Volume I: Elementary Theory and Methods[END_REF]. A MPP provides a naturally sparse representation of configurations of geometrical objects of interest which are inside a considered continuous space. They allow the physical process to be better approximated by a disconnection from the numerical model induced by a regular sampling grid (pixels in image or voxels in volume).

Let K ⊂ R d (d = 2 or 3 in this paper) be a compact subset with 0 < ν(K) < ∞, ν a Borel measure. A realization of a PP on K is a finite configuration of points {k i ∈ K, i = 1, ..., N } as k i = k j for i = j. To form more complex objects, a geometrical mark can be attached to each point, like a radius for example. In that case, a configuration of circles whose points {k i , i = 1, ..., N } are the centers is obtained. Let (M, M, P M ) be the probability space which describes the marks (M ⊂ R + for a radius). A finite random configuration of marked points (or objects) is a sample of a MPP if only the position process of objects is a PP. Based on this definition, volume or image geometrical features are viewed as a set of objects identified jointly by their positions in the image and their geometrical characteristics (the marks).

The simplest MPP is the homogeneous Poisson MPP. That is to say, the number of points in each A ⊂ K has a Poisson distribution with mean ν(A) = βµ(A), with µ the Lebesgue measure, and the points representing the positions of the objects are uniformly and independently distributed in K. The marks are chosen independently in M for each point. For a more complete presentation of MPP, the reader is referred to [START_REF] Stoica | A gibbs point process for road extraction from remotely sensed images[END_REF][START_REF] Chatelain | A bayesian marked point process for object detection. application to muse hyperspectral data[END_REF].

Let us define a configuration of marked points by y = {(k 1 , m 1 ), ..., (k N , m n(y) )}, where k i ∈ K and m i ∈ M , i = 1, . . . , n(y) denote the object positions and the object marks respectively. n(y) is the number of marked points in the configuration y. A configuration of a MPP is classically viewed as a sample issued from an unnormalized probability density f (as the normalizing constant is intractable) which is a Gibbs distribution:

f (y|θ) ∝ exp(-U (y|θ)) (1)
with θ a set of parameters. The energy U (y|θ) allows the properties of configuration y to be modeled and is composed of the sum of two terms: 1) a data-driven energy or external energy denoted U ext (y|θ ext ) that reflects the adequacy between configurations of objects and the observed data and 2) an internal energy denoted U int (y|θ int ) that reflects an a priori on such configurations. This leads to the following expression:

U (y|θ) = U ext (y|θ ext ) + U int (y|θ int ) (2)
and θ = θ ext ∪ θ int . Thus, for a given value of θ, the most likely configurations correspond to the ones that have the global minimum of the total energy:

y = argmin y U (y|θ) (3)
The computation of the global minimum of the energy is classically performed by a simulated annealing which is a stochastic method of optimization (see [START_REF] Geman | Stochastic relaxation:gibbs distributions and the bayesian restoration of images[END_REF][START_REF] Azencott | Simulated annealing: Parallelization techniques[END_REF] and section 2. 3).

An appropriate definition of the data energy allows marked points (objects) that are consistent with respect to a given observation to be obtained. For Tomo-PIV, the process must converge to an appropriate configuration of 2D or 3D objects with respect to the projection images. This section presents our main theoretical contribution in object reconstruction based on MPP: a datadriven energy derived from a minimum mean square error (MMSE) criterion. In this paper, the data-driven energy will be used to search 2D objects in projection images (K ⊂ R 2 ) or to reconstruct 3D objects in a volume (K ⊂ R 3 ), from projection images. In both cases, it is computed thanks to the intensity of 2D objects in projection images. Its definition takes into account a possible interaction of objects when associated geometrical marks overlap (see Fig. 2). Let us notice that previously defined MPP in a context of image processing never introduce an interaction in the data-driven term: the interactions between objects are defined inside the internal term [START_REF] Stoica | A gibbs point process for road extraction from remotely sensed images[END_REF][START_REF] Stoica | Detection of cosmic filaments using the candy model[END_REF][START_REF] Chatelain | A bayesian marked point process for object detection. application to muse hyperspectral data[END_REF][START_REF] Alata | Grouping/degrouping point process, a point process driven by geometrical and topological properties of a partition in regions[END_REF]. In [START_REF] Mallet | A Marked Point Process for Modeling Lidar Waveforms[END_REF], the MMSE is already used to define the data-driven energy of the proposed MPP but the contribution of each marked point and their interaction terms are not developed.

Energy of marked point process

Data-driven energy

The mean square error (MSE) of a configuration y knowing an observation set o = {o s } s∈S can be written as follows:

M SE(y |o ) = s∈S (o s -p s ) 2 (4)
S represents the set of sites and p s the reconstructed value in s by the configuration y:

p s = yj ,yj →s p yj →s (5)
Thus, p s is the sum of the contributions of each object y j ∈ y which has a non-zero intensity at site s. This property will be denoted y j → s and p yj →s is the intensity generated by y j at site s (see Sections 1 and 3.1 for discretization procedure and objects used in Tomo-PIV). Equation 5 also implies that objects have a finite extension (associated to the chosen geometrical marks) and they do not therefore provide intensity all over the sites.

In view of these considerations, a neighboring relationship between two objects can be defined as follows:

y j o ∼ y k if ∃s ∈ S s.t. y j → s and y k → s (6)
The development of M SE(y |o ) (Eq. 4) and the deletion of constant terms (computation details can be found in Appendix A) allow the expression of a datadriven energy in the form of a sum of two energies U ext,1 (y|θ ext ) + U ext,2 (y|θ ext ) to be obtained:

U ext,1 (y|θ ext ) = yj ∈y φ ext,1 (y j ) (7)
where

φ ext,1 (y j ) = s∈S,yj →s p yj →s p yj →s -2o s (8)
and

U ext,2 (y|θ ext ) = yj o ∼y k ,j<k φ ext,2 (y j , y k ) (9)
where

φ ext,2 (y j , y k ) = s∈S,yj →s & y k →s 2p yj →s p y k →s (10)
By considering Equations ( 8) and ( 10), the datadriven energy can be written:

U ext (y) = yj ∈y φ ext,1 (y j ) + yj ,y k ,j<k yj o ∼y k φ ext,2 (y j , y k ) (11)
The first-order term, φ ext,1 , contains a correlation operator between the values generated by an object and the values at its corresponding locations in the observed data. If φ ext,1 (y j ) < 0, the object y j is attractive and its presence allows the global energy of the system to be minimized. Conversely, if φ ext,1 (y j ) > 0, the object y j will be penalized. Then, U ext,1 allows the creation of objects in regions with high intensities in observed images. Note that since some objects are attractive, it is necessary to have an exclusionary term to avoid the accumulation of objects in the final configuration. It is the role of the second-order term U ext,2 . Thus, this energy is based on the neighborhood relationship o ∼ (see Eq. 6) that occurs between two objects when their projections overlap (Fig. 2): function φ ext,2 acts as a correlation operator on the intersection of objects. It is important to note that its computation depends on the spatial extension of the projections of objects and their non-zero overlaps.

To sum up, for a given object y j , φ ext,1 (y j ) is high when the values in the observed images are low, causing an increase in the U ext value. Conversely, when the values in the observed images are high, φ ext,1 (y j ) is negative which favors the corresponding position of the associated object. φ ext,2 (y j , y k ) penalizes two objects which intersect in an image. The data-driven energy will then have low values for a set of correctly positioned objects which do not have too many similar projections. By minimization, it seems to lead to a desired configuration of objects. Now, the second term of the energy is presented: the internal energy.

Internal energy

In the proposed model, the internal energy function is divided into a sum of two terms and can be written as follows:

U int (y|θ int ) = U e (y|β) + U s (y|γ a ) (12)
The first term U e (y|β) = -n(y) log(β) is an energy associated with the intensity of the process in terms of the number of objects n(y) in a configuration. It is defined by the β parameter, which represents the mean number of points by unity of volume in an homogeneous Poisson PP (see section 2.1).

The second term U s (y|γ a ) = -n a (y) log(γ a ) allows a Strauss point process which belongs to the family of Markov point processes to be defined [START_REF] Van-Lieshout | Markov point processes and their applications[END_REF][START_REF] Alata | Grouping/degrouping point process, a point process driven by geometrical and topological properties of a partition in regions[END_REF]. When 0 ≤ γ a < 1, this component penalizes the aggregation of objects and to go further in its definition, we need to briefly describe the objects used for tomo-PIV reconstruction. The objects are Gaussian intensity distribution with limited spatial extension defined by a fixed energy percentage (see Eq. 27 & 28 and section 3.1.1 for more details): in 2D, the spatial shape of an object is a circle; in 3D, it is a sphere. In both cases, the spatial extension of the Gaussian intensity distribution is characterized by a radius. If the geometrical mark does not have such a shape, some other choices could be made for the definition of U s (see [START_REF] Stoica | A gibbs point process for road extraction from remotely sensed images[END_REF] for example).

n a (y) represents the number of neighbor relationships between objects in the following sense:

y i S ∼ y j if ||k i -k j || 2 ≤ r i +
r j where r i and r j are the radiuses of y i and y j respectively. The value of the hyper-parameter γ a ∈ [0, 1] controls the outcome of the potential function. If γ a = 1, the process defined by U int behaves as a homogeneous Poisson PP with intensity β. If γ a ∈]0, 1[, pairs of objects with a distance less than r i + r j are penalized. If γ a = 0, the process forbids that two objects exist within a distance r i + r j . The process is then said to be hard core.

The proposed model is then parametrized by θ = θ ext ∪ θ int with θ int = {β, γ a }. θ ext is mainly defined by the way objects are projected and discretized on images (see Section 1 & 3.1).

OVRMPP -Object Volume Reconstruction based on Marked Point Process

This section describes the MPP dedicated to the reconstruction of 3D objects in Tomo-PIV which represents a particular implementation of the previously proposed MPP. Our aim is to reconstruct a volume of 3D objects based on a set of projection images. Unlike classical tomographic reconstruction methods, the proposed tomographic reconstruction approach, called OVRMPP, allows 3D objects that belong to a continuous 3D space to be obtained (see Fig. 2). In view of the elements recalled in section 2.1, points are center positions of 3D objects and marks provide the set of features that characterize each object in the volume: a radius of a sphere and a central intensity (see section 3.1.1). For a configuration, y = {(k 1 , m 1 ), ..., (k N , m n(y) )}, with

k i ∈ K and m i = {r i , E i } ∈ M , i = 1, . . . , n(y), we have -K ⊂ R 3 . -M ⊂ R +2
M SE(y |{I i , i = 1, ..., N } ) = N i=1 s∈Si (o i,s -p i,s ) 2 (13)
and equations 8, 9 and 10

φ ext,1 (y j ) = N i=1 s∈Si,yj →s p yj →s p yj →s -2o i,s (14)
U ext,2 (y|θ ext ) = N i=1 yj o i ∼y k ,j<k φ i ext,2 (y j , y k ) ( 15
)
where

φ i ext,2 (y j , y k ) = s∈Si 2p yj →s p y k →s (16)
and y j oi ∼ y k means that 2D projections of 3D objects y j and y k intersect in projection image I i . The projection values p yj →s of object y j are computed using the central energy E j and equations 25 and 26. The radius that defines the finite extension of y j is given in equation 27.

In section 2.4, the 2D case used to obtain a configuration of 2D objects for each projection image will quickly be described.

Simulation of MPP and simulated annealing

Metropolis-Hastings-Green algorithm

PP and MPP are classically simulated using the Reversible Jump Markov Chain Monte Carlo method (RJM-CMC) exploiting a Metropolis-Hastings-Green (MHG) dynamics (or algorithm) [START_REF] Green | Reversible jump markov chain monte carlo computation and bayesian model determination[END_REF][START_REF] Stoica | A gibbs point process for road extraction from remotely sensed images[END_REF][START_REF] Perrin | 2d and 3d vegetation resource parameters assessment using marked point processes[END_REF][START_REF] Weina | Marked point processes for crowd counting[END_REF]. This method allows the simulation of a process with varying sampling spaces. In our case, these sampling spaces are associated with configurations with different numbers of objects. The basic moves of RJMCMC for MPP are the birth and death moves [START_REF] Preston | Spatial birth-and-death processes[END_REF][START_REF] Baddeley | Stochastic geometry models in high-level vision[END_REF][START_REF] Lieshout | Markov point processes and their applications[END_REF]: during theses moves, one object is proposed to be either added to or removed from the current configuration.

In addition to these basic moves, supplementary moves can be used: a configuration can be changed by moving positions of objects or by modifying the mark. To obtain a better position, an object inside the configuration is randomly chosen and new positions are randomly sampled using a uniform probability inside an area around the original position. In case of the use of multiple propositions, the proposed position is the one that gives the minimal value of data-driven energies.

Then, the proposed state for a population y containing d-dimensional objects, is y = y ∪ {y p } for a birth move, where y p is a randomly proposed object using a continuous uniform distribution in K and M , if no other prior distribution exists. For the death move, the proposed state is y = y \ {y i }, where y is uniformly chosen inside y. For the translation move, the proposed state can then be written as y = y \ {y i } ∪ {y p } where y i , i ∈ {1, . . . , n(y)} is a chosen object inside y and k p is randomly chosen near k i .

The different propositions are accepted with probability min {1, τ i }, i = B, D or T R (for Birth, Death and Translation), with τ i the MHG acceptance ratio (see Eq. 17, 18 and 19):

if τ i ≥ 1, the proposition is accepted; otherwise, u ∼ U [0,1] is sampled; if u ≤ τ i the proposition is accepted otherwise it is rejected.
Let us denote f b the probability of selecting the birth move and f d the probability of choosing a death move. In the case of a birth move, the position of the new object y is chosen randomly following a uniform distribution on the volume K. Therefore, the MHG ratio can be written as follows when the mark is fixed:

τ B = f d f b f (y |θ) f (y|θ) ν(K) n(y) + 1 ( 17
)
and when the mark is not fixed and a uniform law is assumed, τ B is multiplied by (r max -r min )(I max -I min ) [START_REF] Stoica | A gibbs point process for road extraction from remotely sensed images[END_REF].

In the case of a death move, the point to be removed is chosen with a uniform probability on the existing points in the configuration. The following MHG ratio is obtained when the mark is fixed:

τ D = f b f d f (y |θ) f (y|θ) n(y) ν(K) (18)
and when the mark is not fixed and a uniform law is assumed, τ D is divided by (r max -r min )(I max -I min ).

In the case of a translation move, The following MHG ratio is obtained:

τ T R = f (y |θ) f (y|θ) (19)
(see Sec. 2.1 for ν(K) and f (y|θ) is defined in Eq.( 1)).

The different propositions are accepted with probability min {1, τ i }, i = B, D or T R. In each case, the ratio f (y |θ) f (y|θ) = exp {-∆U }, with ∆U = U (y |θ) -U (y|θ), only depends on the terms associated to the added, deleted or translated object. For example, in the case of a birth:

∆U B = φ d,
where n yp is the number of neighbors of y p in y following the S ∼ neighbor relationship (see Sec. 2.2.2).

MAP estimation based on simulated annealing

As previously stated, the MHG algorithm allows configurations of points to be sampled depending on the defined probability density functions. In order to obtain a configuration that minimizes the energy (see Eq.

3), we need to perform a stochastic method of optimization like simulated annealing. This technique is based on the simulation of a non-homogeneous Markov Chain (see [START_REF] Stoica | Detection of cosmic filaments using the candy model[END_REF][START_REF] Alata | Grouping/degrouping point process, a point process driven by geometrical and topological properties of a partition in regions[END_REF] for examples with MPP). Simulated annealing needs the definition of a new energy which contains a term T which can be viewed as a temperature (see Eq. 11 and 12):

U T (y |θ ) = U e (y |β ) + U ext (y |θ ext ) + U s (y |γ a ) T (21)
and is done by simulating the modified process from a high temperature T 0 to a low temperature T f . U e is not divided by T in order to keep constant the prior about the mean number of objects by volume unity during the iterations. Indeed, if U e is divided by T , it means that β becomes β 1 T and when T goes down to zero, β 1 T will decrease to zero if β < 1 or increase to the infinite if β > 1 which is not suitable.

One of the main interests of simulated annealing is that it is supposed to converge whatever the initial configuration which can be the empty configuration for example. But, in our case, we need the convergence of configurations with more than 50,000 3D particles: the dimension of the space to explore is huge. By decreasing the temperature slowly (the optimal decreasing scheme is also intractable), the algorithm will need a huge computation time. Our objective is therefore to obtain an "accurate" configuration in a reasonable time. To this purpose, we propose an initialization method in order to obtain a first configuration of 3D objects.

Let us note that the simulated annealing provides us with a Maximum A Posteriori (MAP) estimation of the 3D objects: the external energy can be used to define the likelihood function and the internal energy the prior about the solution.

The complete algorithm (called IOD-OVRMPP) contains the following two steps (see Fig. 3): 1) the IOD method (see next section) that provides a set of 2D detected objects for each projection, i = 1, . . . , N and the initial 3D object set computed from the 2D object set; 2) OVRMPP method that provides the final 3D object configuration based on a simulated annealing procedure.

IOD -Initialization method based on Object Detection

In the tomographic reconstruction field where data size is very high, the initialization process is a very important step that can speed up the convergence of the reconstruction methods and enhance the reconstruction quality. To compute a first estimation of the reconstructed volume, we propose an "Object Oriented" method called IOD (Initialization by Object Detection) which was inspired by a 3D position reconstruction technique named "triangulation" often used in the computer vision field. This initialization procedure will be used before the 3D object reconstruction process, whose aim is to obtain a 3D object set with high position accuracy. The operating mode of this 3D object volume initialization method is applied to N acquired projection images and can be carried out in three main steps (Fig. 3): 1) detection of 2D objects in each projection image; 2) identification of epipolar 2D objects; 3) reconstruction of 3D objects from 2D epipolar objects marks.

Fig. 3 Main steps of the proposed method.

1 st step: 2D object detection

To detect 2D objects in the projection images, a MPP for each projection image is used. In view of the elements recalled in section 2.1, points are now center positions of 2D objects and marks provide the set of features that characterize each object in a projection image: a radius of a circle and a central intensity (see section 3.1.1). For a configuration of 2D objects of the ith MPP, i ∈ {1, . . . , N }, y = {(k 1 , m 1 ), ..., (k N , m n(y) )}, with k j ∈ K and m j = {r j , I j } ∈ M , j = 1, . . . , n(y), we have: For the equations of the MPP, see section 2.2 and equation 25 for the computation of the data-driven energy.

-K ⊂ R 2 . -M ⊂ R +2 and
To limit the number of iterations and speed up the convergence of the simulated annealing, we develop a simple initialization procedure, called IRW (Iterative Random Walk), in order to propose a first distribution of 2D objects. To this purpose, one iteration of this procedure is performed in three steps: a simple peak detection is applied to the image which gives a set of detected 2D positions ; each position of this set will be transformed into a 2D object by randomly moving the corresponding position in a short range, in order to enhance its data-driven energy; only a finite number of positions is tested, without the optimization procedure here. A residual image which will be used in the next iteration, is computed between the image used for peak detection and the one generated from the 2D object set. In section 3.2, the number of iterations is fixed equal to two. After this initialization procedure, a simulated annealing optimizes the 2D objects distribution (see Sec.

2.3.2).

2 nd step: 2D epipolar object identification

Once the detection of 2D objects is completed, a search algorithm to identify 2D epipolar objects is applied. To do that, a geometric correspondence between a point X 0 in physical space and a point x i 0 in the projection image I i is determined by the camera optical arrangement. It is described by a projection function ψ i defined as follows:

ψ i (X 0 ) = x i 0 (22)
It is considered that the projection of a point in the volume onto a given image is unique. A back-projection function can also be defined; however, for a point x i 0 in the image I i , there is not one unique corresponding point X 0 in the volume but a set of points. This set of points, in the volume whose projection is x i 0 , is a curve called the line-of-sight (LOS) i (x i 0 ). All the linesof-sight can be parameterized using the Z variable [START_REF] Thomas | Optimization of the volume reconstruction for classical tomo-piv algorithms (mart, bimart and smart): synthetic and experimental studies[END_REF]. For a given Z, the back-projection function is defined as follows:

ψ -1 i (x i 0 , Z) = (X, Y, Z) (23)
In practice, the projection and back-projection functions are obtained through a rigorous calibration procedure [START_REF] Wieneke | Volume self-calibration for 3D particle image velocimetry[END_REF]. The camera model used is a pinhole model but other models could be used.

In a perfect case, a set of N 2D objects, which belong to their respective projections, are considered epipolar neighbors if their lines-of-sight intersect in one point. In reality, there is always some noise in the 2D object positions and the lines-of-sight do not intersect. In order to identify epipolar neighbors, the volume where the intersections can be found is decomposed into 3D boxes, placed on a list. For each box, and for each projection, a list of 2D objects whose lines-of-sight cross the box is identified. When in a given box, the list contains one and only one 2D object the procedure is stopped for this box. If one of the 2D object lists is empty, the box is removed from the list. Otherwise, the box is cut into eight boxes which are added to the list. The algorithm goes on like this until it reaches the end of the list. At the end, the box list contains only boxes with epipolar neighbors candidates. Their "intersection" can be defined as the point which minimizes the average distance with all the lines-of-sight. Some of them are false candidates. Two criteria are used to eliminate them:

their "intersection" must be inside the box the intersection error (minimal average distance between the lines-of-sight) must be lower than a given threshold.

If the noise is too strong, wrong epipolar candidates are identified and there is no way of separating false from true candidates. In this case, it is probably not possible to get a good volume reconstruction. And the higher the particle density, the higher the noise influence is. A box list example is shown in figure 4.

3 rd step: 3D objects initialization

The position of the 3D objects is obtained as the "intersection" of the lines-of-sight of the epipolar neighbors.

The mark is chosen such as the average projection size on all images is the same as the average 2D objects size. The link between the 2D and 3D object marks is detailed in section 3.1.1.

In the next section, the application of the proposed "object oriented" approach to Tomo-PIV field is developed.

Application to fluid mechanics field

In this section, the formalism presented above is applied in the field of fluid mechanics to reconstruct the volumes of 3D particles from a set of projection images. For the sake of brevity, only the experiments made on synthetic data are presented here. We will present some additional results we have obtained studying a turbulent channel flow [START_REF] Earl | Volumetric measurements by tomographic piv of an open channel flow behind a turbulent grid[END_REF] in a further paper.

Tomo-PIV technique

As mentioned in the introduction, during the last decade, a lot of developments has been proposed for improving the reconstruction accuracy and for accelerating the calculation. Part of these results has been obtained in the collaborative European framework AFDAR (Advanced Flow Diagnostics for Aeronautical Research [START_REF] Afdar | Advanced flow diagnostics for aeronautical research[END_REF]). The objective of AFDAR was to focus on new threedimensional methods based on Particle Image Velocimetry (PIV) for measurement in more industrial applications.

The size of the measurement volume within the laser sheet is fully flexible and can be adapted to measure small volumes of a few mm 3 as well as large volumes if the light power is sufficient. However, the depth of focus remains a limiting factor. Flow-tracing particles have a real diameter of the order of 5 -80 µm and a spherical shape. Pixel intensities in images constitute spots that represent 2D particles. 2D structures within images are assumed coming from Mie scattering phenomenon [START_REF] Chýlek | Absorption and scattering of light by small particles: the interference structure[END_REF]. In fact, the diffraction of the laser beams from spherical 3D particles enables its energy to be spread to get 2D projected particles using different viewing angles.

Data description and objects of interest

The projection images (Fig. 5) used to reconstruct volumes are characterized by their high resolution and high 2D particle density (overlapping of particles phenomenon). The density of particles in the images is measured in particles per pixel (ppp). Intensities of pixels in the images are presented in gray level with particles uniformly distributed on a black noisy background (Fig. 5). The image of a particle in the flow is a spot. The intensity in the image is maximum where the probability of finding the particle position projection is the highest. The spot is characteristic of the uncertainty about the particle position. In real experiments, the distance between the cameras and the measurement volume is quite large, which allows the light beams coming onto the sensors of the cameras to be considered quasi-parallel. Under these assumptions, the model can be simplified and the uncertainty about the position of the 2D particles can be considered isotropic. It is also assumed that it follows a normal distribution.

Volumes to be reconstructed are characterized by their large size and their high density of 3D particles. The density of particles in the discrete volume is measured in particles per voxel (ppv).

The main philosophy of the proposed approach is to reconstruct only objects of interest hence the name "object-oriented" method.

In section 2, we presented MPP based on 2D or 3D objects, each one defined by a mark containing an intensity and a radius. This mark is used to compute data-driven and internal energies (see Sections 2.2.1 and 2.2.2). The radius represents the uncertainty about the position of the object. In the next sections, we provide more details about the marks for 2D and 3D objects and the relationship between them.

3D particle

As explained previously, rather than reconstruct the "real" particles, their position and the uncertainty about their position, which is the mark of the objects, are obtained. This approach is linked to our objective to be close to classical methods in order to cross-correlate volumes of particles in order to obtain the most probable displacement.

A 3D particle1 ξ can be identified by a 3D continuous position X 0 = (X 0 , Y 0 , Z 0 )

T ∈ K and characterized by its marks m ξ containing its maximum value E 0 and its radius R. These two parameters define completely the uncertainty about the 3D particle position. It is assumed to be an isotropic multivariate Gaussian law and can be computed using the following equation at

X = (X, Y, Z) T ∈ [-R, R] 3 : E ξ (X) = E 0 exp - (X -X 0 ) T (X -X 0 ) 2σ 2 ξ ( 24
)
where σ ξ is the standard deviation of the 3D Gaussian distribution. Outside a sphere of radius R, E ξ (X) is supposed equal to zero. R is defined with respect to σ ξ and E 0 (see Eq. 27) through a coverage factor. R defines the confidence radius. The projection of the 3D particle uncertainty is supposed to be conservative, in the sense that the uncertainty about the particle projection position is the projection of the 3D probability distribution on the image. So, for small uncertainties, the 3D isotropic Gaussian distribution modeling the presence probability density is projected onto a 2D isotropic Gaussian distribution, that represents the uncertainty about the projection position.

2D particle

Let us define ζ a 2D object obtained by the projection of the 3D object ξ. ζ has the position x 0 = (x 0 , y 0 ) T which belongs to a subset of R 2 . We assume that a 2D particle 2 ζ is characterized by a confidence circle of radius r associated to an isotropic Gaussian intensity distribution with intensity I 0 : m ζ = (r, I 0 ). The coverage factor is the same as for the 3D particles. Therefore, the intensity distribution can be computed using the 2D isotropic Gaussian function defined at x = (x, y) T such as (x -x 0 ) 2 + (y -y 0 ) 2 ≤ r 2 as follows:

I ζ (x) = I 0 exp - (x -x 0 ) 2 + (y -y 0 ) 2 2σ ζ 2 (25)
where I 0 is the center intensity and σ ζ the standard deviation of the 2D Gaussian probability distribution.

Outside the confidence sphere, I ζ (x) is neglected. r is defined with respect to σ ζ and I 0 (see Eq. 28) with the coverage factor.

Projection/Back-projection process

Data driven energy used in our reconstruction approach is computed between images and projection of 3D particles. To do this, functions that allow 2D particles to be obtained from a 3D particle and vice versa are needed. These functions are called the projection and backprojection process. Projection of a 3D particle on N projection images gives N 2D epipolar particles. Backprojection of N epipolar 2D particles gives one 3D particle.

Let ξ be a 3D particle with an isotropic Gaussian probability distribution, identified by its 3D position X 0 and characterized by a central intensity E 0 and a confidence radius R. Let ζ be a 2D particle with an isotropic Gaussian probability distribution, identified by its 2D position x 0 and characterized by a central intensity I 0 and a confidence radius r. Projection and back-projection processes are performed by finding the correspondence between 3D and 2D marks.To get a 2D particle from a 3D particle in the projection image I i , its 2D position x i 0 should firstly be computed using equation [START_REF] Descamps | Automatic flamingo detection using a multiple birth and death process[END_REF].

Once the central position x i 0 of a 2D particle is obtained, its radius and its central probability density can be searched for by evaluating the projection of a 3D particle with isotropic Gaussian distribution profile on 2 If the mathematical notations of section 2 which are given for a generic population of objects are considered, ζ belongs to a population of 2D objects (see Section 2.4.1) and ζ = y i for example. Then, we have the following correspondence: k i = x 0 , r i = r and E i = I 0 a plane.The relation between 2D and 3D central probability density is given by the following equation:

I 0 = E 0 √ 2π σ ξ (26)
σ ζ is proportional to σ ξ . The proportionality constant depends only on the camera model. The standard deviation σ ξ is used to compute a discrete map of a 3D particle to discretize the volume of particles after reconstruction. The discrete maps are generated so that the probability of the particle of being inside a sphere of radius R reaches a confidence level δ. The radius is given by equation [START_REF] Gamboa | Méthode du maximum d'entropie sur la moyenne et applications[END_REF] where k 3D (δ) is the 3D coverage factor:

R = k 3D (δ)σ ξ (27)
For example, for δ = 90% we have k 3D (δ) = 2.55. After the completion of the projection process of a 3D particle, a discrete map is also computed for each 2D particle using the same confidence level δ as for an isotropic Gaussian distribution. Discrete maps related to 2D particles are used to compute the data-driven energy (see Eq. 8 and 10). To get the discrete map, the standard deviation σ ζ of the Gaussian should firstly be computed. In the same way as above, σ ζ is given by:

r = k 2D (δ)σ ζ (28)
where k 2D (δ) is the 2D coverage factor.

To study the performance of the proposed approach, tests were conducted on simulated data (results on real data are given in [START_REF] Ben-Salah | Reconstruction de volumes de particules par processus ponctuels marqués[END_REF]). Our approach was adapted to the measurement of flows and called IOD-PVRMPP (Iterative Object Detection -Particle Volume Reconstruction based on Marked Point Process) since it was used to reconstruct the population of 3D particles. A comparison between IOD, IOD-PVRMPP and MinLOS-MART which is a reference algebraic reconstruction method in Tomo-PIV will be presented in the following sections. The comparison is performed using several evaluation criteria.

Results on synthetic data

A set of synthetic volumes with associated projected images are generated using an image generator developed in C++ using the SLIP library [START_REF] Tremblais | Slip : Simple Library for Image Processing (version 1.0)[END_REF] and inspired by Lecordier and Westerweel [START_REF] Lecordier | The EUROPIV Synthetic Image Generator (S.I.G.)[END_REF]. The synthetic volume size is 500×500×150. Four projections of size 500×500, with image seeding density varying from 0.0004 to 0.2 ppp (Fig. 6), are then computed knowing that 0.05 ppp is the required seeding density level for Tomo-PIV field [START_REF] Elsinga | Tomographic particle image velocimetry[END_REF]. 2D and 3D particles are characterized by marks as described in section 3.1 and [START_REF] Ben-Salah | 3d particle volume tomographic reconstruction based on marked point process: Application to tomo-piv in fluid mechanics[END_REF][START_REF] Ben-Salah | Particle volume reconstruction based on a marked point process and application to tomo-piv[END_REF]. The number of 3D particles N p in synthetic volumes, generated for each seeding density level, vary from 100 up to 50000 particles. For this experimental study, the intensity center and the radius of all particles are fixed: the diameter size is 5 voxels for 3D particles and 3 pixels for 2D particles; the distribution of intensities around all center positions of particles is modeled as a multivariate isotropic Gaussian density. The calibration model is a pinhole model, without distortions. The acquisition system is simulated with 4 cameras (N = 4): two cameras on one side with a 30 • viewing angle; two other cameras on the other side of the volume, on the same plane (Fig. 7). For IOD and PVRMPP, a set of parameters is chosen for each seeding density in order to get a good reconstruction of 3D particle volumes (see tables 1 and 2). These parameter choices were adopted after several simulations, trying to find the best trade-off between reconstruction quality, ghost particle rate and computation time. In tables 1 and 2, f tr corresponds to the probability of selecting the translation move. Birth and death moves are chosen with the probability 1-f tr , and, in a second random selection, with the probabilities f b and f d respectively (f d = 1 -f b ). γ a is fixed to 0.05 (see Section 2.2.3). The simulated annealing algorithm is configured with an initial temperature T 0 and with a fixed final temperature T f equal to 0.02. A classical cooling scheme is chosen: T t = T 0 q t with t the current iteration, q = T f T0 1 N it the parameter of the cooling scheme and N it the number of iterations. Each iteration consists in one proposed "Birth"/"Death"/"Translation" move of MHG algorithm, which may be accepted or not (see Section 2.3). Table 2 Sets of parameters used for the PVRMPP reconstruction method applied to synthetic data.

To evaluate the performances of IOD-PVRMPP method, some quantitative measurements are computed and compared to those obtained with MinLOS-MART [START_REF] Thomas | Optimization of the volume reconstruction for classical tomo-piv algorithms (mart, bimart and smart): synthetic and experimental studies[END_REF], with a fixed number of iterations equal to 5, which can be considered as a reference algorithm within our application. The results of IOD method are also given to show the improvement obtained between initial and final populations of points. All measures are computed by taking 10 reconstructed samples (volumes) for the IOD-PVRMPP method given the random nature of the algorithm.

First, two examples of evolutions of the cumulative mean of the number of 3D objects (N p ) in the population during PVRMPP are presented, with an initial population obtained with IOD (Fig. 8). These curves are truncated showing about the first half of the iterations. Two main types of evolution are observed: for ppp lower than 0.05, N p globally increases before reaching a value close to the true number of particles. For ppp higher than 0.1, IOD provides a population with a lot of generated ghost 3D objects. Then, many 3D objects are deleted at high temperature (during the first itera- tions). After this phase, N p increases and converges to the true number of particles.

To evaluate the reconstruction quality, the crosscorrelations (Q v ) between the reconstructed volumes and the reference volumes are computed:

Q v = Nv i=0 V r (i) -V r V ref (i) -V ref (V r -V r ) 2 (V ref -V ref ) 2 (29)
with N v the number of voxels, V r the reconstructed volume and V ref the reference volume. Figure 9 shows that IOD-PVRMPP provides better results than MinLOS-MART for seeding densities up to 0.1 ppp. The reconstruction quality of the IOD-PVRMPP method decreases when the overlapping between 2D particles in images increases. To evaluate the accuracy of the IOD-PVRMPP method, error positions (E v ) between 3D particle positions in reconstructed volumes and the reference volumes are computed. Results show a very high level of accuracy with error average varying between 0.01 to 0.043 voxels for seeding densities up to 0.05 ppp (Fig. 10). Curves in (Fig. 11) show very low rates of detected ghost particles with IOD-PVRMPP which vary between 0 and 0.098% for seeding density up to 0.05 ppp against 2.92 to 47.74% with MinLOS-MART. The results obtained above can be visually observed in a slice of the reconstructed volumes. In order to view 3D particles inside volumes and visually compare the results, Y -axis slices of volumes reconstructed with 0.05 ppp image density using MinLOS-MART and IOD-PVRMPP methods are extracted. A slice of the reconstructed volume using MinLOS-MART method shows several particle position errors like ghost 3D particles (red circles) and false-negative 3D particles (blue circles) as opposed to the IOD-PVRMPP reconstructed slice that is very close to the reference slice (Fig. 13). The features of 3D particles located in the reconstructed volumes have a direct influence on the quality of the reconstruction and the cross-correlation to generate 3D velocity field. In order to see the influence of the reconstruction quality on the quality of velocity fields without considering correlation algorithm errors, the measurement criterion Q c named correlation quality is used [START_REF] Thomas | Optimization of the volume reconstruction for classical tomo-piv algorithms (mart, bimart and smart): synthetic and experimental studies[END_REF]. This criterion is defined as follows:

Q c = N vf i=0 u ref (i) -u a (i) 2 N vf i u r (i) -u a (i) 2 (30)
where N vf is the number of positions where the velocity field is computed, u ref the velocity field generated from the reference volumes (synthetic), u a the analytical velocity field, and u r the velocity field computed from the reconstructed volumes.

To do this, data with 0.05 ppp level density is used. Particles belonging to each pair of successive volumes follow a simulated sinusoidal movement. The

(U, V, W ) = x d -x ∆t , y d -y ∆t , z d -z ∆t (31)
with ∆t the temporal step fixed at one voxel.

The IOD-PVRMPP method provides a velocity field of quality equal to Q c = 0.99 compared to Q c = 0.89 for the MinLOS-MART method. To generate velocity fields, a multi-pass cross-correlation algorithm developed in the context of the French VIVE3D project using the SLIP library is used [START_REF] Tremblais | Slip : Simple Library for Image Processing (version 1.0)[END_REF]. The cross-correlation results were obtained using search windows with dimension 32×32×32 voxels and an overlap of 75%. Vector fields were validated by a 5×5×5 voxel local median filter size. False vectors were replaced by the median value.

In order to assess the results obtained by calculating the quality correlation criteria Q c , L 2 distance map D between the velocity field and the reference one were computed for MinLOS-MART and IOD-PVRMPP. The results are shown on Fig. 14. They confirm the numerical result obtained with the Q c criterion. Representation of particles via marks takes up very little space compared to representation inside a 3D array volume. Indeed, if, as in our implantation only six parameter values are associated to a particle, the memory size of a set of marks characterizing a 3D particle is equal to 496 bytes. Consequently populations of particles varying from 100 to 50000 elements used in our experiments need only 0.05 to 24.80 MB storage. In comparison, for the MinLOS-MART method, the required memory is linked to the dimensions 500×500×150 voxels of the discrete volume which corresponds to 286 MB storage. Gain in memory storage for IOD-PVRMPP increases as the size of the volume increases.

Conclusion

In this paper, we have presented a new approach called IOD-OVRMPP for the reconstruction of a population of 3D objects, including an initialization procedure, using the Marked Point Process framework. Rather than a "voxel-oriented" approach, our work uses an "objectoriented" approach. It allows the problem to be solved in a parsimonious way and provides a MAP estimation: the reconstruction of a 3D object set is obtained by minimizing an energy function thanks to a simulated annealing based on an MHG algorithm. To speed up the convergence of the simulated annealing, an initialization method has been developed. It provides the initial distribution of 3D objects based on the detection of 2D objects located in projection images. The IOD-OVRMPP approach works in a continuous space where 2D and 3D objects are identified by real positions. This mathematical formalism is very suitable for the reconstruction of sparse solutions since it reconstructs only objects, unlike voxel-oriented methods of tomographic reconstruction which use all the voxels to finally provide a discrete volume.

The proposed approach is applied to the Tomo-PIV field to reconstruct volumes of 3D particles. The method is validated on synthetic data. Results show the relevance of the proposed data-driven energy for the reconstruction of a volume of 3D particles. The method is compared to an algebraic reconstruction method called MinLOS-MART. The IOD-PVRMPP method improves the quality of reconstruction even for high particle densities (up to and even more than 0.05 ppp). The position error and ghost particle rate are greatly reduced with IOD-PVRMPP: 3D position error is equal to 0.043 voxels, compared to 0.25 voxels with MinLOS-MART for 0.05 ppp particle density. For the same particle density, the number of ghost particles is approximately 0.1%, against 47.7% with MinLOS-MART. IOD-PVRMPP also provides a better memory usage: the amount of RAM used to run the method is less compared to classical volumetric memory usage. Future work will include the optimization of the IOD-PVRMPP algorithm to speed up the reconstruction as well as a full comparison with the Shake the Box method [START_REF] Schanz | ???Shake The Box???: A highly efficient and accurate Tomographic Particle Tracking Velocimetry (TOMO-PTV) method using prediction of particle positions[END_REF] or the SMTE method for time resolved acquisitions [START_REF] Lynch | An efficient and accurate approach to MTE-MART for time-resolved tomographic PIV[END_REF]. Then the IOD-PVRMPP method will be improved in order to handle particle blurring due to a small depth of field and astigmatism due to imaging optics. This approach will also be tested in other domains such as digital holography for example, to see its capacity to reconstruct 3D object volumes using different data.

Appendix A Data Driven Energy calculation details

If the square of equation ( 4) is developed, the following result is obtained:

M SE(y |o ) = s∈S o 2 s + p 2 s -2o s p s (32)
The first term of the addition in equation ( 32) is constant. The data-driven energy defining the likelihood between the projection of a configuration y and the observed data o = {o s } {s∈S} can be written:

U ext (y) = s∈S p 2 s -2o s p s (33)
U ext defines the quality of a configuration compared to the data: the closer the projection of a population of objects is to the reference image, the lower its energy value. We show below that this energy can be written as a sum of first-order and second-order neighborhood energy terms. Using equation [START_REF] Azencott | Simulated annealing: Parallelization techniques[END_REF], equation ( 33) can be written: [START_REF] Lecordier | The EUROPIV Synthetic Image Generator (S.I.G.)[END_REF]

U ext (y) =

Fig. 1

 1 Fig. 1 Tomo-PIV principle.

Fig. 2

 2 Fig.2Overlapping of two 2D objects issued from the projection of two separate 3D objects, located on the same line of sight.

Fig. 4

 4 Fig.4Example of a box list. The colored lines represent the lines-of-sight, the black spheres, the theoretical particles and the orange circles, the lines-of-sight "intersection" from noisy data.

Fig. 5

 5 Fig. 5 Real image of 2D particles with 1200×1600 pixel size and an estimated seeding density of 0.035 ppp.

Fig. 6

 6 Fig. 6 Bands extracted from simulated images with different concentrations of particles. The density of particles of each band is indicated in particles per pixel (ppp).

Fig. 7

 7 Fig. 7 Camera positions for the simulated cases.

Fig. 8

 8 Fig. 8 Cumulative means of the number of particles inside the volume during two PVRMPP runs. top: ppp = 0.05; true number = 12500; initial number after IOD = 7925; estimated number = 12491. bottom: ppp = 0.1; true number = 25000; initial number after IOD = 25704; estimated number = 24448.

Fig. 9

 9 Fig. 9 Evolution of the volume reconstruction quality Qv as a function of the density of particles.

Fig. 10

 10 Fig. 10 Evolution of the mean error position inside the volume as a function of the density of particles.

Fig. 11

 11 Fig. 11 Evolution of percentage of number of ghost particles as a function of the density of particlesSome samples of projections after volume reconstruction are presented in figure[START_REF] Benson | Modified simultaneous iterative reconstruction technique for faster parallel computation[END_REF] in order to provide some qualitative results.

Fig. 12

 12 Fig. 12 Comparison of the projections (region with size 64×64) of a reconstructed volume of particles onto a camera for (a) MinLOS-MART, (b) Reference (c) IOD-PVRMPP for a 0.05 ppp seeding density level. The images are originally encoded with 16 bits unsigned integers and differences between the images (b) and (c) are difficult to observe.

Fig. 13 Y

 13 Fig. 13 Y -axis slices issued from volumes reconstructed with MinLOS-MART and IOD-PVRMPP methods using 0.05 ppp image density. Green circles represent true particles, red circles represent ghost particles and blue circles represent nonreconstructed particles (false-negative particles). The magenta color is derived from the superposition of red and blue circles.

  3D coordinates (x d , y d , z d ) of a particle after displacement depend on its initial coordinates (x, y, z) and are calculated using the equations x d = A sin(ω z) + C cos(ω y), y d = B sin(ω x) + A cos(ω z) et z d = C sin(ω y) + B cos(ω x) with A = √ 3, B = √ 2, C = 1 and ω = 1/20. The coordinates (U, V, W ) at each position of the velocity field depend on the coordinates of the particles after the motion and are calculated as follows:

DFig. 14 L 2

 142 Fig. 14 L 2 Distance map D (mm/s) on XY plane, to assess IOD-PVRMPP and MinLOS-MART methods for 0.05 ppp level density particles.

→s p yj →s p yj →s -2 o s + s∈S 2 yj ,y k ,j<k yj →s,y k →s p yj →s p y k →s ( 36 )

 236 k ,j<k yj →s,y k →s p yj →s p y k →s U ext (y) = yj ∈y s∈S yj →s p yj →s p yj →s -2 o s f irst-order term + yj ,y k ,j<k yj o ∼y k s∈S yj →s & y k →s 2p yj →s p y k →s second-order term

  and minimum and maximum values for intensities (I min and I max ) and radiuses (r min and r max ) need to be fixed in order to completely define M .-The observation data is the set of the N projection images {I i , i = 1, . . . , N } with I i = {o i,s } s∈Si , S i the set of pixels of I i . -The overlapping of 3D objects must be strongly penalized. Then γ a (see Section 2.2.2) needs to be taken close to 0 (see Section 3.2 for the values used).

	Equation 4 becomes:

  1 (y p ) +

	ny p a
	yj ∈y,yp

o ∼yj φ d,2 (y p , y j ) -log βγ

  minimum and maximum values for intensities and radiuses need to be fixed in order to completely define M . -The observation data is a projection image: I i = {o i,s } s∈Si , S i the set of pixels of I i . -The overlapping of 2D objects often occurs. Then γ a (see Section 2.2.2) needs to be taken close to 1 (see Section 3.2 and 3.2 for values used).

Table 1

 1 Sets of parameters used for the IOD method applied to synthetic data.

If the mathematical notations of section

which are given for a generic population of objects are considered, ξ belongs to a population of

3D objects (see Section 2.2.3) and ξ = y i for example. Then, we have the following correspondence:k i = X 0 , r i = R and E i = E 0

Acknowledgements The current work has been conducted as part of the AFDAR project, Advanced Flow Diagnostics for Aeronautical research, funded by the European Commission program FP7, Grant No. 265695 and also the FEDER project No. 34754.


Gaël Privault

Guénolé Huitric

Marius Hervé

Elzbieta Trzop

Julien Tranchant

Benoit Corraze

Zohra Khaldi

Laurent Cario

Etienne Janod

Jean-Claude Ameline

Nicolas Godin

Roman Bertoni
email: roman.bertoni@univ-rennes1.fr

Ultrafast photo-induced dynamics of V2O3 thin films under hydrostatic pressure

Keywords: Ultrafast Spectroscopy, Photo-Induced Phenomena, Thermo-Elastic effects

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Nowadays, a new field of research is related to the ability of controlling the functionality of materials by using intense external field. Indeed, the next challenge is not only to observe and to understand physical processes leading to the emergence of interesting properties in materials (magnetism, conductivity, etc) but to directly control the desired functions in the fastest and most efficient way. From that perspective, the current availability of ultrafast intense laser pulses has created a pivotal change in the last decades. It is indeed possible to act on the different degrees of freedom (spins, orbital, lattice and electrons) using electromagnetic waves covering a broad range of energy. Moreover, thanks to femtosecond laser pulses, the light-matter interaction may occur at the relevant time scale and induce nonthermal distributions. Such concept is extensively developed in the field of photo-induced phase transition where an intense light pulse is used to generate new phases such as photo-induced superconductivity [1,2], spin transition states [3,4] or even "hidden phase" [5]. Those ultrashort light pulses are also relevant to investigate fundamental couplings in condensed matter such as electron-electron [6][7][8] or electron-phonon [9][10][11][12]. When dealing with Photo-Induced Phase Transition (PIPT), the temperature is very often the only external thermodynamical parameter used to explore the phase diagram of materials. This severely restricts the attainable area in the phase diagram and prevents thorough investigation of photo-induced dynamics near transition lines. From that perspective, the complementary use of hydrostatic pressure opens interesting new perspectives related to understanding of non equilibrium dynamics near critical end points in phase diagram. V 2 O 3 is a system often referred to as prototypical Mott material displaying a metal to insulator transition (MIT) (Figure 1) [13][START_REF] Hansmann | Mott-Hubbard transition in V2O3 revisited[END_REF][START_REF] Mcwhan | [END_REF]. Under ambient conditions with a small amount of chromium used for chemical substitution, the material is in the paramagnetic insulating phase (PI) with an electronic gap of hundreds of meV. By applying pressure, the system enters in the paramagnetic correlated metallic phase (PM) as the electronic bandwidth is modulated and the electronic gap collapse. In this system, an antiferromagnetic phase (AFI) is present at low temperature with a transition temperature around 150 K at ambient pressure. According to the phase diagram of V 2 O 3 , one does not expect to cross any transition lines by applying pressure on pure V 2 O 3 . The system is indeed in the PM phase under ambient conditions and the application of hydrostatic pressure only drives the system further in the PM phase. Nevertheless, the application of hydrostatic pressure may be reflected in the out-of-equilibrium response of the material after photo-excitation. By using ultrafast pump-probe measurements, it is possible to follow the dynamics induced by femtosecond optical excitation on the relevant time scale. From that aspect, broadband optical spectroscopy is a useful tool that may provide valuable information regarding the complex out-of-equilibrium dynamics generated by femtosecond optical excitation [16].

Results

Mirror Window Substrate Film

R Tr

Fig. 2 Sketch of the experiment. Pump and probe beams are represented in dark red and khaki respectively. Both transmitted (Tr) and reflected (R) beam can be collected for measurement.

Here, we report the transient evolution of optical properties of a thin film of V 2 O 3 after ultrafast photo-excitation by femtosecond laser pulses. The experimental observables are made of broadband optical spectra (450-725 nm) measured in two configurations: reflectivity and transmission (Figure 2). Both pump and probes beams arrive onto the film in quasi-normal incidence (less than 4 • ) after passing via a hollowed parabolic mirror. The thin film is sandwidched by the optical saphirre window and the substrate (0.5 mm thick). The sapphire window and the thin film are not perfectly parallel so reflected beams can be spatially separated. The pump laser energy is set to 1.55 eV (800 nm) with an estimated penetration depth of about ∼ 90 nm [17] while the film thickness is approximatively 190 ∓ 10 nm. The first step of this study is to perform ultrafast optical spectroscopy in the pressure cell at ambient pressure (1 bar) while the material is already surrounded by helium gas. The transient transmission data obtained at ambient pressure is shown in figure 3. Photoinduced transient changes over the visible spectra last more than two hundreds of picoseconds. By slicing the spectrum at given wavelengths, one can observe more clearly the ongoing photoinduced dynamics (Figure 4). First, an instantaneous peak occurs around time zero where pump and probe overlap in the time domain. This ultrashort response relaxes very quickly. Such a feature can be ascribed to a mixture of ultrafast electronic response [18] and non-linear optical effects [19] induced by the intense electric field of the pump in the transparent sapphire windows. Consequently, a thorough investigation of those early dynamics often related to electron-electron dynamics is difficult in this experimental configuration. Afterwards, the transmission increases up to ∼ 30 ps for all probed wavelengths. After passing this maximum, the signal undergoes an oscillating behaviour before reaching a plateau around ∼ 100 ps. By looking at different wavelengths (Figure 4), one can observe that only the amplitude of signals changes but the overall photo-induced dynamics remain similar. This data set present a chirp correction that sets a proper timing for all wavelengths as explained in the section 4. A clear oscillating component modulating the full spectrum is observed and last for more than one picosecond. It is better seen when looking at a single wavelength dynamic around 650 nm (Figure 6). The frequency analysis of these oscillations indicates a single component centred around 7.3 THz (inset of Figure 6).

Ambient pressure measurements

Tr

The transient reflectivity data are shown in figure 7. A small ultrashort peak lasting less than 1 ps around time zero is present for some wavelengths with the same physical origin than the one observed in transmission geometry. After, the different dynamics are observed with an initial negative signal for all wavelengths during the first tens of picoseconds. Later, a slower oscillating signal appears with an average value that differs depending on the wavelength (Figure 8). Nevertheless, similar to the transmission case, this slow oscillating behaviour occuring after 40 ps seems to be rather unaffected by the probing wavelength. This measurement in reflection geometry is very sensitive to acoustic effects and propagation of strain wave [20,21]. It is sometimes referred to as time-resolved Brillouin scattering (TRBS) [22].

High pressure measurements

In order to look for the impact of hydrostatic pressure onto the photo-induced dynamics, ultrafast measurements are performed under a hydrostatic pressure of 6 kbar. All other experimental parameters remain similar compared to the ambient pressure case. for ambient pressure. One can see the dynamics in both geometries are very similar to those observed at ambient pressure. The transmission data reveals a transient increase of signal up to ∼ 20 ps followed by a small oscillating pattern. transient signal reach a plateau after 80 ps. The transient reflectivity data set presents a negative signal lasting tens of picoseconds followed by an oscillating pattern around a mean value. In general, both signals follow the same trends pointing out that the very same photo-induced effects are generated in the high pressure case.

Discussion

The material being in the high temperature PM phase, the interaction between the optical pump pulse and V 2 O 3 results in photo-excitation of electrons generating an excited electronic gas within the pump pulse duration (60 fs). First, the electronic subsystem is going to thermalize on sub-picosecond time-scale to reach a Fermi-Dirac statistical distribution [7,18,23]. Then, it interacts with the lattice through electron-phonon coupling in order to attain a thermal distribution in both the electronic and phononic subsystems. This exchange energy process should occurs on few to tens picosecond time-scale depending on the materials [9,10,12,18] and is often modeled by a two-or multi-temperature model [24]. Nevertheless, it does not imply that the lattice has truly reached a thermal distribution [25]. The resulting effect is an ultrafast lattice heating with an in depth profile matching the optical penetration depth of the pump pulse. This process results in the generation of strain waves along the film thickness that affect the refractive index of the material [START_REF] Huitric | Faraday Discussions[END_REF]. The resulting dynamics already reported in several materials belongs to the class of photo-induced thermo-elastic effects that can be probed through optical spectroscopy [21,[START_REF] Huitric | Faraday Discussions[END_REF][START_REF] Levchuk | [END_REF][28][29] or X-ray diffraction [30][31][32]. Despite being already few decades old and starting with the seminal work of Thomsen [20], this field of research is still vivid thanks to the development of new laser and X-ray sources and apparition of new techniques. The penetration depth profile of the 26.5×10 -12 = 7.15/km.s -1 in agreement with value of 7.3/km.s -1 reported in literature for ambient pressure [START_REF] Yelon | Solid State Communications[END_REF]. The slight shift towards earlier value (∼ 25 ps) the high pressure case support a modification of the speed of sound in the material with an estimated value of ∼ 7.6/km.s -1 . This small but detectable shift is better seen in figure 11 where both transmission signals are directly normalized and compared. Also, the absence of clear wavelength dependence in the observed signals (Figure 4 and8) independently of the applied pressure rules out the effect of time-resolved Brillouin scattering as been responsible for the oscillating behaviour observed up to 100 ps. Much likely, this effect is related to the travelling back and forth of strain waves inside the thin film that may be affected by the acoustic impedance with the substrate. Those oscillations are call acoustic eigen-modes in case of free standing films [START_REF] Huitric | Faraday Discussions[END_REF][START_REF] Levchuk | [END_REF]. The plateau lasting for hundreds of picoseconds afterwards is related to the deformation of the thin film induced by laser heating that will relax on longer time-scale via heat diffusion and energy exchange with the substrate. The reflectivity measurements show different behaviours than the ones observed in transmission despite originating from the same photo-induced effects. The first negative signal directly comes from interferences phenomena between the pulse reflected from the surface and the one partially reflected by the strain wave propagating into the material. This effect, in the case of semi-infinite material, is refer as time domain Brillouin scattering [20,[START_REF] Parpiiev | [END_REF]. It also explains the strong wavelength dependency expected for TRBS [22] and observed within the first 20 ps (Figures 8 and10). Indeed, the first positive extremun is shifted from 19 to 14 ps when changing wavelengths from 700 to 550 nm at ambient pressure. On contrary, the latter oscillating component lasting up to hundreds of picoseconds does not exhibit a clear temporal shift while changing the probing wavelength. This apparent absence of temporal shift is an indication that the observed dynamic does not correspond to TRBS. It seems more close to the observation of acoustic eigen-modes of the film whose frequencies should by independent of the probing wavelength [START_REF] Levchuk | [END_REF]. When comparing transient reflectivity observed under ambient and high pressure (Figure 12), similar conclusions than the ones made for transmission can be drawn. A blue shift of the oscillating feature is observed in the high pressure case suggesting a more rigid material with increased speed of sound. Those results point out that even a moderate pressure of 0.6 GPa is enough to modulate the thermo-elastic photo-response of a thin film of V 2 O 3 . When looking at the sub-picosecond dynamics, oscillatory components are always present irrespective of the experimental geometry or the hydrostatic pressure (Figure 13). Both curves have been sligtly shifted for clarity. It reveals that the excited electronic gas preferentially coupled to the zone center fully symmetric A 1g optical mode to release excess energy. The measured value of 7.3 THz (Figure 6) matches the ones reported in previous experiments [18,[START_REF] Huitric | Faraday Discussions[END_REF]35]. Under the investigated pressure range, no modification of the phonon frequency are observed as shown in figure 13, neither detectable variation of damping time. It implies that this phonon potential remains almost unperturbed within our sensitivity range up to 0.6 GPa which is expected for a rigid inorganic material such as V 2 O 3 .

Methods

Sample preparation

A thin film of pure V 2 O 3 is deposited on a c-cut sapphire substrate by sputtering of a vanadium target in Ar discharge. Afterwards, the thin film is annealed at 530 • C in a controlled atmosphere yielding in a single phased V 2 O 3 layer. The thickness of 190 ± 10 nm is estimated via scanning electronic microscopy.

Optical measurements

All the results of this study are based on the classical pump-probe method. An intense ultra-short light pulse at a given photon energy is used to photo-excite the system (pump) while a second one of lower intensity (probe) is used to monitor the evolution of the system over time. The probing beam is made of a supercontiuum generated in a sapphire plate [START_REF] Alfano | The supercontinuum laser source: fundamentals with updated references[END_REF]. A typical spectrum is shown in the upper panel of figure 14. The spectral cut in the UV side comes both from the coupling between the probe pulse and the optical fibre used to bring the supercontinuum into the spectrometer and the drop of reflectivity of silver mirrors. The sharp edge in the IR side comes from a low pass filter used to remove noise from the pump pulse centered at 800 nm. We can also estimate the sensitivity of the setup by performing the analysis implemented for time-resolved data in the absence of pump pulse. The result is shown in the lower panel of figure 14, displaying a global sensitivity of 2.10 -4 (0.2 mOD) from 475 nm (2.6 eV) to 725 nm (1.7 eV). The pump photon energy is centred at 1.55 eV (800 nm) with an estimated duration of 60 fs. The typical size (FHWM) of the pump pulse of quasi-gaussian profile is around 160 µm leading to an average fluence of 50 µJ.mm -2 . This fluence is used for all experiments reported in this study. Each transient time-resolved spectra consists in an average of multiple single shots of excited or ground states. As the pump is chopped at half the laser frequency (500 Hz), alternative pump and unpumped spectra are measured by the camera. The spectra are sorted thanks to a small leak of the pump pulse onto the CDD camera that is used as a marker for excited spectrum. After, for each delay, 5000 shots of excited and unexcited spectrum are average and used to extract transients transmission and reflectivity. We define the transient transmission as the following: Where T R (of f, λ) corresponds to the transmission of the unperturbed material and T R (t, λ) to the one at a given time delay (t) after the photo-excitation. In the same way, the transient reflectivity in function of time is defined as followed:

∆T R (t, λ) = T R (t, λ) -T R (of f, λ) T R (of f, λ) (1)
∆R(t, λ) = R(t, λ) -R(of f, λ) R(of f, λ) . (2)
Where R(t, λ) is the reflectivity of excited states and R(of f, λ) the reflectivity of the unperturbed ground state. Those two quantities are the ones used to monitor the evolution of the material after ultrafast photo-excitation.

For very short delays, we may correct the data set for the chirp experienced by the probe beam. This effect arises from the ground velocity dispersion experienced by the probe beam and implies that each wavelength creates its own timing referential. Nevertheless, if the spectral resolution is sufficient, the main limitation for the experimental time resolution comes from the pump pulse duration [START_REF] Polli | [END_REF]. This is illustrated by the figure 15 showing chirp uncorrected and corrected data set where the A 1g optical phonon of V 2 O 3 is clearly visible in both cases despite the 1.25 ps duration of supercontinuum pulse in the probed range. It emphasises that all photo-induced informations are already embedded in data sets without chirp correction. The complete high pressure setup is built on a multi-stages amplification process. First, a conventional gas compressor is used to pressurize the entire setup up to 1000 bars. Secondly, the primary stage is isolated and a piston-based compression is used to boost up the cell to 7 kbars. The high-pressure gas is delivered to the cell via a metallic capillary. A schematic transverse cut of the cell is drawn in figure 16. The pressure is monitored at the entrance of the capillary, in the primary stage. It allows for direct monitoring of the He gas pressure in-operando with precision of few bars. This is a substantial advantage over conventional diamond-like anvil cell relying on indirect measurement such as ruby fluorescence. The accuracy of the current apparatus is ± 10 bars (at 7 kbars) and ± 1 K in the operating range as the cell fits into a home-made cryostat for cooling. It possesses two sapphire windows (Figure 16) used for both reflectivity and transmission optical measurements. The use of sapphire optical windows ensures good transmission in a spectral range spanning from 5 eV to 0.25 eV (250 to 5000 nm). The sample area, depicted in red, has a volume of 0.4 cm 3 allowing for the investigation of thin films as large as 5 × 5 mm 2 . In addition, this gas cell apparatus, controlled by electro-valve, allows ramping and relaxing pressure with arbitrarily slow rate, as low as few bars per minute. This setup is especially relevant for the study of materials undergoing first order transition accompanied by significant volume changes. Indeed, a very slow variation is needed to accommodate the strain experienced by the material during the transition.

Conclusion

In conclusion, we demonstrate the ability to perform ultrafast optical spectroscopy in a high pressure He gas cell. The current setup is able to reach time resolution as good as ∼100 fs while conserving a signal to noise ratio of 10 -4 . Using V 2 O 3 as a benchmark material, the obtained results show a clear effect on the transient modulation of optical properties of a thin film of V 2 O 3 mainly mediated by acoustic effects. Indeed, despite being a rigid inorganic material, clear modifications of dynamics induced by thermo-elastic effects are observed under the moderate pressure of 0.6 GPa. Nevertheless, we do not observed obvious modifications of the coherent optical A 1g mode suggesting only a small modification of the lattice potential under the investigated pressure. The possibility of very fine tuning of pressure could reflect on very accurate control of the photo-acoustic response of material. This capability is particularly relevant in the emerging field of photo-induced strain driven phenomena where the generation and control of strain wave could be used to induce phase transition [32]. Furthermore, this new experimental platform should be relevant for the study of soft correlated molecular materials where complete phase transition occurs within few kbars range allowing a large exploration of phase diagram [38].

Acknowledgments. Authors acknowledge Agence Nationale de la Recherche (ANR) for funding undergrant, ANR-21-CE30-0011-01 CRITICLAS. This work was carried out in the frame of the International Laboratory, DYNACOM IRL University of Tokyo -CNRS -UR1.

Fig. 1

 1 Fig. 1 Schematic phase diagram of V 2 O 3 . Red cross denotes the ambient thermodynamical conditions. Blue arrow illustates the expected displacement over the phase diagram.

Tr

Fig. 3 Fig. 4

 34 Fig. 3 Transient transmission of V 2 O 3 thin film after photo-excitation with 800 nm pump under ambient pressure

Fig. 5

 5 Fig. 5 Transient transmission of V 2 O 3 thin film after 800 nm excitation on sub-picosecond timescale under ambient pressure

Fig. 6

 6 Fig. 6 Transient transmission of V 2 O 3 on sub-picosecond at 650 nm. inset) FFT of the oscillating signal

RFig. 7 Fig. 8

 78 Fig. 7 Transient reflection of V 2 O 3 thin film after photoexcitation with 800 nm pump

Fig. 9 Fig. 10

 910 Fig.9Transient transmission of V 2 O 3 thin film after 800 nm excitation at selected wavelengths under 6 kbar

Fig. 11

 11 Fig. 11 Comparison of the transient transmission at 650 nm between 1 bar and 6 kbar

Fig. 12

 12 Fig. 12 Comparison of transient reflectivity at 650 nm between 1 bar and 6 kbar

Fig. 13

 13 Fig.13 Transient oscillating dynamics related to A 1g phonon at 1 bar and 6 kbar.

Fig. 15

 15 Fig. 15 Sub-picosecond data set without (top pannel) and with (lower pannel) chirp correction

Fig. 16

 16 Fig. 16 Transverse cut of the pressure gas cell with relevant dimensions. The two sapphhire windows are shown in blue. The shaded red space denotes the sample area.

Declarations

• The authors declare no conflict of interest.