Ningning Song

Carole Deumié

Anabela Da
email: anabela.dasilva@fresnel.fr

Considering sources and detectors distributions for quantitative photoacoustic tomography

Keywords:

Considering sources and detectors distributions for quantitative photoacoustic tomography References and linksNoninvasive laser-induced photoacoustic tomography for structural and functional in vivo imaging of the brain

Introduction

Photoacoustic tomography (PAT) is a multiwave imaging modality that potentially combines the high optical contrast brought by the optics and high spatial resolution by ultrasound probing. The interests in using this technique for biomedical sensing are growing since 1990s, and its applications are in breast cancer detection [START_REF] Ermilov | Laser optoacoustic imaging system for detection of breast cancer[END_REF][START_REF] Heijblom | Visualizing breast cancer using the Twente photoacoustic mammoscope: What do we learn from twelve new patient measurements?[END_REF], small animal molecular imaging [START_REF] Razansky | Multispectral photoacoustic imaging of fluorochromes in small animals[END_REF][START_REF] Li | Photoacoustic tomography and sensing in biomedicine[END_REF], subcutaneous structures imaging [START_REF] Zhang | In vivo imaging of subcutaneous structures using functional photoacoustic microscopy[END_REF], functional imaging [START_REF] Wang | Noninvasive laser-induced photoacoustic tomography for structural and functional in vivo imaging of the brain[END_REF], and other applications.

The principle consists in illuminating a sample with a local absorbing optical contrast (vasculature, presence of a tumor for example). The absorbed energy converts into heat and the temperature of the object increases, causing thermal expansion that generates an acoustic pressure wave. The energy deposited is proportional to the fluence and the local absorption coefficient of the tissues. Hence, the primary aim of PAT is to retrieve the spatial 3D distribution of the absorption coefficient of the probed tissue. This absorption coefficient, proportional to the concentration of the different chromophores composing the tissue, is an intrinsic marker of the physiology and the metabolism. In the red and near infrared, PAT has been applied to mapping hemoglobin concentrations that can reveal the presence of a cancer, through oxy-and deoxy-hemoglobin concentrations which are markers of the oxygen consumption in the tissues [START_REF] Zhang | Imaging of hemoglobin oxygen saturation variations in single vessels in vivo using photoacoustic microscopy[END_REF]. Glucose monitoring with PAT has also been reported [START_REF] Kottmann | Glucose sensing in human epidermis using mid-infrared photoacoustic detection[END_REF]. However, usually the technique is applied to probing tissues at small depths (small animals, subcutaneous examination), for which the fluence can be modelled as a Beer-Lambert's law by considering light scattering is not perturbing too much the quantification. For deep tissues imaging, because of their high level of scattering, diffusion of light has to be taken into account in the evaluation of the fluence that depends non-linearly on both the absorption and diffusion coefficients. Isolating the absorption contribution from the energy deposited is the challenge in quantitative PAT (QPAT), as the fluence does depend itself on the optical parameters (absorption and diffusion coefficients) distributions.

In a first approximation, PAT can be limited to a purely acoustic problem: from timeresolved measurements collected at the periphery, acoustic sources are localized through the reconstruction of the initial pressure distribution, from which one can obtain the absorption coefficient distribution by making the assumption of a simple photon propagation model or a uniform distribution of light. The problem is similar to seismology passive sensing: one seeks at localizing the sources that produced the collected acoustic signal. For all these approaches, a wide field illumination is preferable in order to produce a homogenized initial pressure. To go further in the quantification, different strategies have been adopted taking advantage of specific experimental protocols.

There exist basically four approaches in QPAT which consists in recovering:

i. the absorption coefficient directly. Supposing the diffusion coefficient is known, Ripoll et al. [START_REF] Ripoll | Quantitative point source photoacoustic inversion formulas for scattering and absorbing media[END_REF] recovered the absorption map with a point source accounting for instrumental factors such as the source strength, the shape of the optical pulse, and the impulse response and finite size of the transducers; Cox et al. [START_REF] Cox | Two-dimensional quantitative photoacoustic image reconstruction of absorption distributions in scattering media by use of a simple iterative method[END_REF] reconstructed the absorption coefficient assuming the scattering coefficient is known a priori by use of a fixed-point iterative method; Banerjee et al. [START_REF] Banerjee | Quantitative photoacoustic tomography from boundary pressure measurements: noniterative recovery of optical absorption coefficient from the reconstructed absorbed energy map[END_REF] recovered the absorption coefficient directly from boundary pressure measurements. Very recently, Zhou et al. [START_REF] Zhou | Calibration-free absolute quantification of particle concentration by statistical analyses of photoacoustic signals in vivo[END_REF] have introduced a calibration-free method based on the direct measurement of the counts of particles through the statistical analysis of the photoacoustic signal fluctuations.

ii. the absorption coefficient and the fluence. The fluence contribution can be filtered by making use of a contrast agent of known absorption coefficient with photoacoustic signal measurements performed before and after injection [START_REF] Rajian | Quantitative photoacoustic measurement of tissue optical absorption spectrum aided by an optical contrast agent[END_REF]. Daoudi et al. [START_REF] Daoudi | Correcting photoacoustic signals for fluence variations using acousto-optic modulation[END_REF] proposed to combine photoacoustics with acousto-optics allowing also to overcome the contribution of the fluence in the initial pressure measurement. The fluence map can be explicitly measured through Diffuse Optical Tomography [START_REF] Bauer | Quantitative photoacoustic imaging: correcting for heterogeneous light fluence distributions using diffuse optical tomography[END_REF] (DOT). By introducing prior knowledge on the spatial frequency of the two quantities, Rosenthal et al. [START_REF] Rosenthal | Quantitative optoacoustic signal extraction using sparse signal representation[END_REF] extracted both the absorption coefficient (high frequency) and the fluence (low frequency) from the photoacoustic image. As this algorithm is not based on the explicit solution the theoretical light transport equation, it does not require explicit knowledge of the illumination geometry.

iii. the absorption and diffusion coefficients simultaneously by introducing prior knowledge (spectra of known species, diffusion spectral profile) with multiple wavelength illuminations; Bal and Ren [START_REF] Bal | On multi-spectral quantitative photoacoustic tomography in diffusive regime[END_REF] have proven that, when multiple wavelength data are available, the absorption, diffusion and Grüneisen coefficients can be reconstructed simultaneously under additional minor prior assumptions.

iv. the absorption and diffusion coefficients simultaneously by introducing active probing with multiple optical illuminations. In the last decade, a lot of efforts have been produced to improve the quantification through the uncoupled reconstructions of the different physical parameters contained in the initial pressure distribution map, that is: absorption and diffusion coefficients and the Grüneisen coefficient accounting self-consistently on the acoustic properties of the tissue [START_REF] Bal | Inverse diffusion theory of photoacoustics[END_REF]. Tarvainen et al. [START_REF] Tarvainen | Reconstructing absorption and scattering distributions in quantitative photoacoustic tomography[END_REF] reconstructed both the absorption and diffusion coefficients, by considering both the radiative transport and diffusion equations as light transport models; Bal and Ren [START_REF] Bal | Multiple-source quantitative photoacoustic tomography[END_REF] recovered the absorption, diffusion, and Grüneisen coefficient mathematically when the propagation of radiation is modelled by a second-order elliptic equation; Shao et al. [START_REF] Shao | Estimating optical absorption, scattering, and Grueneisen distributions with multiple-illumination photoacoustic tomography[END_REF] estimated the absorption, scattering and Grüneisen distributions from reconstructed initial pressure map with multiple optical sources; In their updated work [START_REF] Shao | Iterative algorithm for multiple illumination photoacoustic tomography (MIPAT) using ultrasound channel data[END_REF], they recovered the absorption and diffusion coefficients directly from measured acoustic pressure. With this last approach, the probing becomes active and similar with what occurs in geophysics when searching for buried objects, petroleum layers, cracks... In the present work, this approach is adopted and tested under different sources and detectors geometries. The purpose is to test its robustness under potential experimental situations with single wavelength illumination, point scanning single element transducer and limited view angle examination. The paper is organized as follows: Section 2 presents the theoretical background, the results are presented and discussed in Section 3, Section 4 is the conclusion.

Theory

Hereafter are recalled the definitions of the physical quantities involved in the photoacoustic phenomenon and the main lines leading to the expression of the solution of the forward model. The inverse problem is also presented. It was formulated according to the method proposed by Shao et al. [START_REF] Shao | Iterative algorithm for multiple illumination photoacoustic tomography (MIPAT) using ultrasound channel data[END_REF] for multiple illuminations and detections situation.

Photoacoustic phenomenon modelling in biological tissues

The amount of heat generated by tissue is proportional to the strength of the radiation. Within the conditions of thermal and stress confinements, light pulse is typically few ns, the heating function can be treated as an instantaneous function proportional to a Dirac () t δ function:

(,) ()(,) ()() () () () . aa H tt t E t µφ µφ δδ =≈ = rr rr r r (1)
Here (,) H t r is the heating function defined as the thermal energy converted at spatial position r and time t by the electromagnetic radiation per unit volume per unit time. ()

φ r is the light fluence [W.cm -2 ],
()

a µ r [cm -1
] is the absorption coefficient. When large tissues are considered, with illumination in the wavelength range (600-1200 nm) for which light scattering mean free path length is much shorter than the absorption mean free path length, the Diffusion Approximation holds [START_REF] Vo-Dinh | Biomedical Photonics Handbook[END_REF] and the spatial distribution of light fluence can be obtained through the resolution of the following diffusion equation:

()() . ( () () ) () a DQ µφ φ -∇ ∇ = rr r r r (2)
() D r [cm] is the diffusion coefficient and () Q r [W.cm -3
] is the illumination light source density. The absorbed energy converts into heat and the temperature of the tissue increases, and thermal expansion takes place, generating an initial acoustic pressure in the medium expressed as 0 (, ) ( , can be formulated linearly with the following matrix form:

( ) TT δ δ =Δ ⇔ = JJ u J H u b (5)
Here model

Δ=pp

, H is the Hessian. The multiplication by T J allows handling an inversion process with smaller dimensions of the matrices involved (reduced to the dimension of the chosen reconstruction mesh). The main problem remains in expressing J . The medium is first meshed into L voxels of volume

l V Δ . Introducing () (,)
1, s S ∈ , becomes: model 1 (,,) ( ) 4 L s dl s dl l l dl s R pt t E V Rv τ τ δ π = Γ =- Δ  ru (6)
and the Jacobian matrix: β ur represents the sensitivity of the energy instantaneously deposited due to the optical perturbation δ u located at positions j r , j r belonging to the reconstruction mesh:

0 model 0, , 0 1 (() ) (, ,) (()
) L s dj s d d l l s j l pt τ τ τ αβ = = =∂ ∂ =  uu Ju r r u u u r (7) , () 4
0 0 0 00 ,0 0 (,) (, ) () () () (() ) () (,) () () sl sl l j a l s aj la j ls j s lj sl al j E ED D φ φδ µ µ µ β φ µ = = =   ∂ +   ∂    ∂∂   ==  ∂∂    ∂    ∂    
φ = =  ∂ =- - - Δ  ∂    ∂  =∇ - ⋅∇ - Δ  ∂   uu uu ru
rru r r ru r r ru rru r r ru r r [START_REF] Ripoll | Quantitative point source photoacoustic inversion formulas for scattering and absorbing media[END_REF] In order to avoid infinite values of 0 G and φ at positions j l = rr , the reconstruction mesh

( [ ] ,1 , j jJ ∈ r , j V Δ
pixel surface) for the values of the optical properties was chosen to be slightly different from the modelling mesh (

[ ] ,1 , l lL ∈ r
).

Computation of the Jacobian matrix

In the following case studies, the simulations were conducted in 2D. The chosen test object (Fig. 1) was basically the same as in [START_REF] Shao | Iterative algorithm for multiple illumination photoacoustic tomography (MIPAT) using ultrasound channel data[END_REF]: a 6 cm × 6 cm square, in the center area of which a 2 cm × 2 cm square was selected as the reconstruction area (Fig. 1, left). The background absorption and diffusion coefficients are respectively (dimensions: 0.3cm × 1.1cm and 0.3cm × 0.3cm, 0.6 cm from the center). The perturbations were chosen relatively small (10%) in order to fulfil the Born approximation. The values correspond to optical properties measured in various situations with DOT (for instance in breast cancer imaging [START_REF] Ntziachristos | Concurrent MRI and diffuse optical tomography of breast after indocyanine green enhancement[END_REF]). If the perturbations are high compared to the background optical properties, non-linear reconstruction algorithms have to be used [START_REF] Shao | Iterative algorithm for multiple illumination photoacoustic tomography (MIPAT) using ultrasound channel data[END_REF]. The sources are located 0.3 cm away from the reconstruction area, both the number and the length of the sources are variable. Point detectors are located 2 cm from this computation area, right at the periphery of the object. ). The FEM mesh typically contains 4300 elements, the modelling mesh is a regular grid of 61 × 61 pixels (pixel size: 0.1 cm), and the reconstruction mesh is 21 × 21 pixels, 1 pixel shifted in horizontal and vertical with the modelling mesh. ,0 (() )

ls j
β ur is then obtained by assembling the sensitivity matrices, and its dimension is () 2 SL J ×× . α is ( )

M TL ××.
The dimension of the Jacobian matrix J is large () 2 TMS J ××× but the Hessian H dimension is reduced 22 J J × . Vector Δ representing the discrepancy between model and measurements has dimension () 1

TMS ××× and b is only 21 J × . The largest matrices J and Δ are computed once and are not involved in the reconstruction process.

Resolution of the inverse problem

Two different reconstruction methods were tested:

-Singular value decomposition (SVD) method for which the Hessian matrix was decomposed into:

T =× × HUSV . S is a diagonal matrix with same dimension as H containing nonnegative singular values i σ ordered in decreasing order, U and V are unitary matrices. The inversion of (5) becomes: -The algebraic reconstruction technique (ART) [START_REF] Kak | Principles of Computerized Tomographic Imaging[END_REF] that searches iteratively for the solution δ u such that:

( ) 2 1
,

nn n n δδ δ - +   =+ -   uu bH u H H , n is a
hyperparameter that can be tuned in order to regularize or to speed up the reconstruction. In all what follows, n was kept constant and equal to 1. Non negative values constraint is also introduced in the algorithm.

These reconstruction methods are tested for different noise levels. For small noise levels, SVD provides sufficiently good results but is known to be unstable and very sensitive to noise, while the ART behaves better, especially when constrained.

Simulation study on sources and detectors geometries

Sets of synthetic data were generated by solving the forward problem and introducing a fraction of random noise with normal distribution. Hereafter are presented the results of simulations performed on the synthetic phantom (Fig. 1) with the above-described method. The purpose of the study is to highlight the effects of different sources and detectors distributions and geometries on the reconstruction of the perturbation of the optical parameters [] T a δδ δ = u D , in terms of accuracy and robustness to noise. To that end, the following aspects were considered: (i) using different reconstruction methods with different noise level: this preliminary step allowed to select the proper reconstruction algorithm; (ii) number of point sources distribution: the point source illumination has been specifically studied in [START_REF] Shao | Iterative algorithm for multiple illumination photoacoustic tomography (MIPAT) using ultrasound channel data[END_REF], here the results are shown as reference; (iii) point detectors distributions: the purpose is to understand which transducers distribution would be the best, for a given number of measurements; (iv) using extended sources instead of the point sources: point sources have been studied but, in practice, extended sources are easier to handle experimentally; (v) limited angle examination: this corresponds to actual experimental situations, here the transmission geometry, where sources and detectors belong to different half-spaces, is studied.

Different reconstruction methods

Four point sources and sixty detectors (15 detectors on each side of the object) are considered for this study. A first point that can be stressed is the poor conditioning of the Hessian: after performing the SVD, and examining the spectrum of the singular values i σ of H, the condition number obtained for this configuration is max( ) min( )

ii cond σσ = = 1.46 × 10 10 . Reconstructions are expected to be unstable and highly sensitive to noise.

The two different reconstruction methods, SVD (with regularization when indicated) and ART (with non-negative constraint, 1000 iterations, and 1, n n λ =∀) were tested under different noise levels (0 up to 10 -2 ). Figure 2 (Left) shows the maps of the reconstructed perturbations on the absorption ( a δ ) and diffusion ( δ D ) coefficients through the different reconstruction methods with increasing noise levels. On Fig. 2 (Right) are represented crossplots extracted from these maps along a horizontal line through the middle of the reconstruction area.

As expected, reconstructions are highly sensitive to noise, a high value of the regularization parameter has to be introduced in the SVD, showing only few elements of H out of noise may be useful for the reconstruction, even at low noise level. With the introduction of non-negative constraint, ART seems to have better potential in the recovery of the optical parameters at higher noise levels. It is interesting to notice that, though ART reconstructs the perturbations with a higher contrast. In all what follows, ART was chosen as the reconstruction algorithm.

Hence, the aim of the following studies is to improve the conditioning of H such that the reconstructions are less sensitive to noise, by optimising the data acquisition geometry.

Sources distributions

In this study, 60 point detectors were distributed around the object and the sources number was varied (2, 4, 8 and 16 sources considered). Examining the condition number of H as a function of the number of sources (Fig. 3), the reconstructions should be indeed more stable as the number of sources increases. Although the inversion is mathematically demonstrated, 2 sources are not enough to obtain uncoupled robust reconstructions.

( ) 2 2
,, /

i j ij ij i j ij QE X X X =-
. The error in the reconstruction values decreases from 1 with 2 sources to 0.4 with 4. An improvement of a factor 10 in the condition number (Fig. 3, Right Top) is obtained by adding only two more sources. However, the improvement in increasing the number of sources might be reduced after 8 sources and the gain in the reconstructions might be small even if the number of measurements (and hence the duration of the examination) is increased. This is all-the-more visible in the cross-plots shown in Fig. 3 Right Bottom: the values reconstructed become closer to the target values, the edges of the objects are sharper and the reconstruction artifacts decrease.

Detectors distributions

Four sources and 60 point detectors are now considered, the number of measurements is now kept constant but the distribution of the detectors is varied. Five situations, schematized in Fig. 4, were examined: a) even and sparse distribution around the object; b) 4-sides examination but short spacing between detectors; c) same situation but with even shorter (compressed) spacing; d) 2-sides examination horizontal detection; e) 2-sides examination vertical detection.

Examining the reconstructions (Fig. 4 Top Left), the quadratic errors (Fig. 4 Bottom Left) and the conditioning of H (Fig. 4 Top Right) obtained for these five situations shows no major improvement in the reconstructed images. Through the values of the condition number and the quadratic errors, one may conclude that distributing the detectors evenly (Fig. 4 situation (a)) may be the most favorable situation. However, reducing the field of view of the detection to the area of interest (Fig. 4 situations (b) and (c)) seems to improve the quantification while not degrading to much the localization. Figure 4 situations (d) and(e) show explicitly the influence of the positioning of the detectors in comparison to the object geometry, especially in the rectangular absorption perturbation: situation (d), when detectors are probing the object through its thinnest width, the emitted acoustic signal is better defined in time and higher in frequency (small objects producing sharper time-resolved acoustic signals), while in situation (e), the resolution is much more degraded in the direction parallel to the detectors (crossplots, Fig. 4 Right Bottom).

Point sources versus wide field illumination

QPAT has been demonstrated by using point sources illumination [START_REF] Shao | Iterative algorithm for multiple illumination photoacoustic tomography (MIPAT) using ultrasound channel data[END_REF]. However, in practice, a wide field illumination is easier and safer to handle experimentally because the irradiation dose can be decreased by distributing the energy over a larger area. For this study, four line sources illuminating the four sides of the reconstruction area were considered (Fig. 5). Their lengths were varied from zero, corresponding to the point source situation, up to the maximum length allowing a full angle probing of the reconstruction area. Increasing the length of the source brings major improvements in the conditioning of H (Fig. 5 Right Top) and, hence, in the reconstructions (Fig. 5 Left Top). These results show that the quantification (see QE Fig. 5 Left Bottom) and localization (see crossplots Fig. 5 Right Bottom) are both definitively improved by using wider field illumination. Hence, the number of measurements required for obtaining accurate reconstructions in terms of localization and quantification can be reduced. When comparing reconstruction results shown in Fig. 3 for situation (c) and those in Fig. 5 for situations (b) to (e), with two times less measurements, the reconstructions obtained are better with 4 line sources illuminations than with 8 point sources (QEs are sensibly the same). However, the medium has still to be probed under multiple points of views: as shown for comparison on Fig. 5 situation (f), using one single wide source does not bring sufficient information to discriminate absorption and diffusion. In this last test case, the absorption coefficient is well reconstructed while the diffusion coefficient is far underestimated, showing that this illumination protocol senses essentially absorption abnormalties (see cross-plots on Fig. 5 Right Bottom).

Example of simulations for experimental situations with limited angle examination

Examining the object under all angles may be experimentally complicated or simply not feasible. Here are reported reconstructions obtained under constrained examination geometry: the medium is probed exclusively under transmission geometry, sources and detectors can rotate and probe the object on its four sides. The purpose of these calculations is to test if satisfying reconstructions can be obtained under less favorable situations corresponding nevertheless to simpler experimental protocols. Figure 6 Top Left shows the five different synthetic experiments considered here. For all of them, the four sides of the object are probed in transmission, sequentially (this corresponds to 4 measurements with an object rotating around its center of gravity with rotation angles 0°, 90°, 180° and 270°), with (a) 1 point source illumination on one side and 15 point detectors evenly distributed on the other side; from this reference situation, the illumination source was enlarged, (b) 1 line source (length: 26 mm) and 15 detectors; then the number of detectors was reduced up to a very limited number (c) same line source and 1 detector only (in this situation, sources and detectors are interchanged compared to situation (a)); (d) same as (c) but with 2 point detectors; (e) same as (c) with 3 point detectors. The condition numbers and quadratic errors obtained for the corresponding reconstructions are also reported Fig. 6.

A first remark is the reconstructions obtained with point or line sources illuminations (Fig. 6 situations (a) and (b)) are improved compared to those for which the detectors are collecting signal from all the sides of the object: when the object is probed in transmission, although the condition numbers are several orders of magnitude higher, still one can get reconstructions of satisfying quality (QE≤0.4) with highly reduced number of measurements (60 measurements only). Figure 6(c) is the oversimplified situation, a single wide source and a single point detector: the absorption coefficient is satisfyingly reconstructed while the diffusion coefficient is not. However, slightly increasing the number of detectors (2 or 3 point detectors per view Fig. 6(d) and 6(e)) allows recovering reconstructions of improved quality.

This result shows the possibility of using simplified experimental PAT setups. Reducing the number of detectors points, up to 3 only here, corresponding to 12 measurements only, of course worsens the quality of the reconstructions but still allows getting better results than with complete angle detection: in comparison with Fig. 5(e), obtained with 240 measurements, the number of measurements is reduced by a factor 20. Additionally, one may think in combining the measurements in order to improve the quality of the reconstructions: the absorption coefficient seems to be fairly well reconstructed with the simplest situation (1 source, 1 detector, Fig. 6 (c)), these results could be used as prior knowledge in a second step reconstruction using 2 or 3 additional measurements (Fig. 6(d) or 6(e)).

Conclusion

In the present work, a multiple illumination photoacoustic tomography algorithm has been implemented according to the method suggested in [START_REF] Shao | Iterative algorithm for multiple illumination photoacoustic tomography (MIPAT) using ultrasound channel data[END_REF]. The method is quantitative in the sense that it allows reconstructing both the absorption and diffusion coefficient while the conventional approaches in PAT, with a single wide field illumination give the access to the initial pressure distribution map. A description with explicit practical implementation details was provided. Different linear reconstruction algorithms were tested: the ART with nonnegative constraints was shown to be less sensitive to noise than SVD or LSQR. Even if not involved in these reconstruction loops, the computation of the Jacobian matrix is still a burden, especially if one seeks to extend the approach to the non-linear case. Gradient-based methods may be used instead because more memory-efficient as they require only gradient information to approximate the Hessian matrix [START_REF] Bal | On multi-spectral quantitative photoacoustic tomography in diffusive regime[END_REF][START_REF] Gao | Bregman methods in quantitative photoacoustic tomography[END_REF][START_REF] Saratoon | A gradient-based method for quantitative photoacoustic tomography using the radiative transfer equation[END_REF].

With this tool, a 2D simulation study on a variety of sources and detectors geometries situations was conducted. The sets of simulations were conducted through the examination of the conditioning of the Hessian matrix. The reconstructions highlighted several important results:

-Optical parameters are reconstructed with improved quality by using wide field illuminations: this study demonstrates that it is preferable to use enlarged sources instead of multiplying the number of point sources illuminations.

-The position of the detectors with respect to the perturbations to be reconstructed is crucial and deserves additional prospections. When comparing the influence of sources and detectors, the role of the illumination and detection schemes are nested and complex. The illumination scheme seems to play a major role in the quantification, while the detection in the localization with high resolved time-offlight measurement. An important result is that probing the sample under transmittance geometry only provides better reconstructions than with full angular coverage detection.

-Different simple experimental situations were tested and good quality reconstructions were obtained with a number of measurements extremely reduced. Experimental validations may support these findings that may influence future experimental setups designs. The present study has highlighted the importance of the data acquisition protocol. A smart combination between different types of measurements will undoubtedly improve the quality of the reconstructions. An example of such a process would be to combine one measurement taken under a full angle illumination, with a field of illumination as wide as possible, and a set of data acquired under different illumination angles. The first measurement would serve in a first reconstruction loop that would provide initialization and constraint values for a second reconstruction loop, performed with the second series of data.

  rr , the model for the considered measurable quantity at time

  evolution transfer matrix propagating the pressure created at initial time at position l r to point detector located at d r ; ,0 (() ) ls j

δG

  is the Kronecker symbol. Within the perturbation approach, knowing the Green's function 0 of the diffusion Eq. (

Fig. 1 .(

 1 Fig. 1. Left: Geometry of the simulated object, light sources and transducers are placed respectively 0.3 cm and 2 cm away from the object. In the reconstruction area are placed the perturbations: two absorbers ( a δµ , positions: (0.6cm;1cm) and (1.6cm;1cm), dimensions: 0.3cm × 1.1cm and 0.3cm × 0.3cm) and one diffuser ( D δ , position: (1cm;1cm), dimensions: 0.5cm × 0.5cm). Right: Different meshes used in the simulations: left, FEM mesh { } Δ r ;

2. 3 . 2 α

 32 Time evolution transfer matrix propagating the pressure , dl τ This matrix represents the acoustic propagation between each pixel of the modelling mesh and the transducer. It carries essentially the information on the distance between detectors and the different pixels of the modelling mesh dl d l R =rr . It is calculated by considering the points located at the center of each pixel of the modelling mesh { } l r . l V Δ is the surface of each pixel, The Grüneisen coefficient is chosen to be a classical one in biological tissues 0steps are considered with a time interval of 0.2 µs. The dimension of , dl τ

1 -

 1 b . A Tikhonov regularization was adopted to avoid instabilities caused by poor conditioning of H : the values of the diagonal of

Fig. 2 .

 2 Fig. 2. Comparison of different reconstruction methods at different noise levels: Left, results of the reconstructions of the perturbations on the absorption a δ and diffusion δ D coefficient maps; Right: cross-plots of the values extracted from a horizontal line in the middle of the vertical axis in the reconstruction areas.

Figure 3

 3 Left Top shows the reconstructions obtained: with 2 sources the diffusion coefficient perturbation is almost not detected, but interestingly, the absorption perturbations are well recovered. The gain obtained in the reconstructions can be quantified through the quadratic error (QE, Fig.3Left Bottom) between the reconstructed values  

Fig. 3 .

 3 Fig. 3. Influence of the number of sources. Left Top: reconstructed perturbation absorption a δ and diffusion δ D coefficient maps with (a) 2, (b) 4, (c) 8 and (d) 16 point sources illuminations. Left Bottom: Quadratic errors on the reconstructions of the absorption (black circles) and diffusion (blue squares) coefficients perturbations as a function of the number of point sources. Right Top: Normalized singular values (Left) and condition number of (H) as a function of the number of sources (Right). Right Bottom: cross-plots of the values extracted from a horizontal line in the middle of the vertical axis in the reconstruction areas.

Fig. 4 .

 4 Fig. 4. Influence of the detectors distributions. Left Top: reconstructed perturbation absorption a δ and diffusion δ D coefficient maps. (a) the detectors are evenly distributed around the object (15 detectors on each side); (b) on each side, the detectors cover a length that is half the length of the object; (c) the length covered by the detectors is half of (b); (d) 30 detectors evenly distributed on each of the two vertical sides of the object; (e) 30 detectors evenly distributed on each of the two horizontal sides of the object. Left Bottom: Quadratic errors on the reconstructions of the absorption (black circles) and diffusion (blue squares) coefficients perturbations for the five different configurations. Right Top: Normalized singular values (Left) and condition number of (H) as a function of the number of sources (Right). Right Bottom: cross-plots of the values extracted from a horizontal line in the middle of the vertical axis in the reconstruction areas.

Fig. 5 .

 5 Fig. 5. Influence of the sources shape. Left Top: reconstructed perturbation absorption a δ and diffusion δ D coefficient maps with (a) 4 points sources; 4 line sources with a length of (b) 0.2 mm; (c) 4 mm; (d) 8 mm; (e) 26 mm; (f) full field of view probing with one single source composed of lines. Left Bottom: Quadratic errors on the reconstructions of the absorption (black circles) and diffusion (blue squares) coefficients perturbations as a function of the length of the four sources ("Integral" represents situation (f), reconstructions obtained with a single measurement but with full field of view source illumination). Right Top: Normalized singular values (Left) and condition number of (H) as a function of the number of sources (Right). Right Bottom: cross-plots of the values extracted from a horizontal line in the middle of the vertical axis in the reconstruction areas.

Fig. 6 .

 6 Fig. 6. Top Left: Schemas of the experimental situations. For each situation, only 4 measurements were performed at 0°, 90°, 180° and 270°, the rotation angles corresponding to rotations around the axis located at point O (cross marker), and perpendicular to the object plan. Schemas correspond to the measurement taken at 0°: (a) 4 points sources, 15 detectors; (b) 4 line sources with length 26 mm, 15 detectors; (c) same line sources, 1 detector; (d) same line sources, 2 detectors; (e) same line sources, 3 detectors; for each situations, four synthetic measurements were considered: these configurations (0°), and object rotated by 90°, 180° and 270°. The condition numbers of the corresponding (H) are reported for each situation. Top Right: reconstructed perturbation absorption a δ and diffusion δ D coefficient maps. Bottom: Quadratic errors (QE) on the reconstructions of the absorption (black circles) and diffusion (blue squares) coefficients perturbations for the five different configurations.

  )

	representing the perturbation to these reference maps. If	0 uu , one can express the δ <<
	measurements as a Taylor series expansion, truncated here at first order (Born
	approximation):									model 00 ,) (, ,) s d pt (, sd p t τ τ δ += ru u ru	+	0 ( ) δ J u u , s d τ
	model ds 0 () s pt (, ,) d τ τ =∂ Ju r u	∂ u	= uu	0	is the Jacobian matrix. Hence, the global forward problem,
	for all considered sources	s		∈	[ ] 1, S	and detectors	[ ] 1, dM ∈	, at any time	[ ] 1, tT τ ∈
			p	ru	=Γ	r	a µφ r ru , Γ being the Grüneisen parameter, dimensionless
	quantity representing the efficiency of conversion of absorbed energy to pressure, and
	u	=	[] T a D . The photoacoustic pressure (, ,) pt ru then propagates at the speed of sound in
	the medium s v as:
					2 2 t  ∂∂ 22 0 (, ,) (, ) ( () ) s vpt p t t δ -∇ =  ∂ ∂  ru ru	(3)
	which solution in free space is well known. For simplicity, to avoid the calculation of time
	derivatives, the measurable quantity considered in this work is the model (, ,) pt ru	defined as a
	quantity directly proportional to the velocity potential [4]:
			model pt v d t p 2 ' 0 (, ,) t s == ' (, , ) t  ru ru	d r	'	0 p 4 π	' (,) ' -ru rr	δ	( t	-	s v -rr	'	)	(4)
	2.2 Inverse problem formulation
	From limited number M of measured pressures, at points i r , due to a limited number S of
	sources, we want to recover the maps of the optical properties	u	=	[	a	] T D . Within a
	perturbation approach, one can always express any optical map distribution as	0 =+ uu	δ	u ,
	where 0 u is the vector representing the maps of the optical properties of a reference medium,
	with known optical properties distribution, and δ u being considered as the vector

Acknowledgments

We gratefully acknowledge China Scholarship Council scholarships for graduate student Ningning Song. This work was partly supported by ANR (AVENTURES-ANR-12-BLAN-BS01-0001-04).


ouverte pluridisciplinaire HAL, est destinée au dépôt

Theoretical study of electronic structure and magnetism in materials for spintronics


Adam Houénou
email: a.houenou@sherpa-eng.com

Philippe Bonnifait

Véronique Cherfaoui

Risk Assessment for Collision Avoidance Systems

published or not. The documents may come   L'archive ouverte pluridisciplinaire

Risk Assessment for Collision Avoidance Systems

Adam Houénou * , Philippe Bonnifait + and Véronique Cherfaoui + Abstract-Collision Avoidance Systems need to perform scene analysis and risk assessment in order to react conveniently. Based on the information provided by the perception system, scene analysis has to predict the evolution of the current driving situation for the near future. Thanks to the predicted trajectories of the relevant traffic participants, the risk of collision on the ego vehicle can be calculated. In many cases, a predicted trajectory is not defined with explicit equations but is given as a set of sampled poses, each one corresponding to a different future time instant. A predicted trajectory being always uncertain, confidence has to be estimated on the so predicted poses. We present a method that propagates the known error covariance matrix of the current pose of the ego vehicle by considering local approximations of the predicted trajectory. This allows to estimate the risk of collision of the ego vehicle with a considered target object. The proposed approach uses a Monte Carlo simulation to approximate the probability that the ego vehicle and the object come into collision at a given future time instant. Each sample time of the whole prediction horizon is considered as a potential collision time so that a curve describing the variation of the risk of collision is obtained. This allows the system to have a better comprehension of the scene and to react proportionally to the threat. The overall approach has been tested with simulated data and the consistency of results is shown.

I. INTRODUCTION

Driver Assistant Systems in general and more specifically Collision Avoidance Systems (CASs) are very active fields of research and development. Recent years have seen many of those systems on the market of series vehicles. A CAS usually performs in low speed scenarios but enhancements are under research in order to cover more scenarios. A CAS needs to predict any eventual collision as soon as possible and must warn the driver or take autonomous actions in order to avoid it or mitigate its consequences. For this purpose, it has to rely first on an accurate perception system that uses sensors, communication devices and data fusion algorithms to detect and track obstacles in the surrounding of the ego vehicle. Then, it is crucial to perform a correct scene analysis which consists in predicting the evolution of the scene. Then, it is possible to assess the risk of collision of the ego vehicle. Depending on the risk assessment result, various criteria can be used to decide the relevant actions that have to be finally executed.

This paper first focuses on the scene analysis part. Based on the estimated current and past states of the scene, the goal here is to understand the current driving situation and to use this result to predict the most probable evolution of the scene. In a previous work [START_REF] Houenou | Vehicle trajectory prediction based on motion model and maneuver recognition[END_REF], we described a trajectory prediction method based on the combination of two predictions made with different approaches. The first prediction uses a deterministic approach. Its assumes a kinematic motion model to predict the trajectory of a vehicle. The second prediction is based on a maneuver recognition algorithm, a trajectory generation process and a cost function used to select the best trajectory. A combination rule allows the final prediction to take advantage of the good accuracy of the motion model approach for short term prediction and of the better accuracy of the maneuver recognition method for longer term prediction. Section II gives a reminder of this prediction method in order to make this paper more intelligible.

Based on this information, the CAS has to evaluate the threat on the ego vehicle in the near future in order to be able to avoid an eventual collision with any potential obstacle. This is the task of the risk assessment module that takes as inputs the predicted trajectory of the ego vehicle and the predicted trajectories of the detected objects. Since these predictions are not deterministic, it is necessary to estimate their uncertainty. A first contribution of this paper is the presentation of a method that propagates the estimated uncertainty of the current pose of a given vehicle along its predicted trajectory that is time sampled and that has no closed-form. This is done by local linearizations of the predicted trajectory and a sequential estimation of the covariance matrices of its consecutive future sample poses. This is described in section III.

Many risk assessment approaches exist in the literature [START_REF] Berthelot | Handling uncertainties in criticality assessment[END_REF], [START_REF] Laugier | Probabilistic analysis of dynamic scenes and collision risks assessment to improve driving safety[END_REF], [START_REF] Lambert | A fast monte carlo algorithm for collision probability estimation[END_REF], [START_REF] Lefèvre | Risk assessment at road intersections: Comparing intention and expectation[END_REF], [START_REF] Polychronopoulos | Dynaminc situation and threat assessment for collision warning systems: the euclide approach[END_REF] but the notion of risk of collision is usually related to the time-to-collision (TTC) or other derived time metrics. The smaller the TTC, the higher the risk. The predicted trajectories being uncertain, the TTC is also uncertain. The proposed approach consists in sampling the prediction time horizon and considering each sample time as a potential TTC. Then, for each sample time and each target object, the probability that the ego vehicle is in collision with the object is calculated. This gives a pseudo-continuous evolution of the risk of collision with each obstacle on the whole prediction horizon which allows the system to define various levels of criticality based on the pairs TTC/probability of collision. The risk assessment method is presented in section IV and its principle is compared to known approaches in the literature. Section V shows the evaluation of the overall risk assessment method in simulated scenarios and the results are discussed. Finally, section VI concludes the paper.

II. TRAJECTORY PREDICTION

From the driver's control commands to the actual trajectory of the vehicle, there are complex interactions including the vehicle's dynamic and the road's structure. Predicting the trajectory of a vehicle (either the ego-vehicle or any detected vehicle) is not a deterministic task because the driver's control commands are quite unpredictable and many parameters are difficult to measure in normal conditions e.g. tire friction on the road. In some works [START_REF] Tamke | A flexible method for criticality assessment in driver assistance systems[END_REF], [START_REF] Polychronopoulos | Dynaminc situation and threat assessment for collision warning systems: the euclide approach[END_REF], [START_REF] Althoff | Lane-based safety assessment of road scenes using inevitable collision states[END_REF], the prediction is made by assuming a kinematic model of the vehicle's motion. Unfortunately, such motion models as constant velocity or constant yaw rate movement are accurate only for a short term, except if the driver's command actually remains unchanged for a long time which is unlikely. During a lane change maneuver for instance, the movement changes a lot and the prediction can quickly be very wrong. In other papers, a Maneuver Recognition Module (MRM) is used to compute a predicted trajectory that takes into account the possible changes in the vehicle's behavior with respect to the recognized maneuver. In most cases [START_REF] Laugier | Probabilistic analysis of dynamic scenes and collision risks assessment to improve driving safety[END_REF], [START_REF] Hülnhagen | Maneuver recognition using probabilistic finite-state machines and fuzzy logic[END_REF], [START_REF] Wiest | Probabilistic trajectory prediction with gaussian mixture models[END_REF], [START_REF] Sung | Trajectory clustering for motion prediction[END_REF], the MRM is based on a set of learning trajectories. The current recorded sequence of the vehicle's trajectory is compared to that database according to various possible techniques and the most likely future trajectory is obtained from the closest trajectories in the database.

In [START_REF] Houenou | Vehicle trajectory prediction based on motion model and maneuver recognition[END_REF], we proposed a new method that mixes both approaches in order to take advantage on the high accuracy of the motion model method for short-term prediction and to respect, for long-term prediction, the possibly changing dynamic of the vehicle during specific maneuvers. Our MRM is not based on a data set; it uses the last measurements of kinematic parameters and the knowledge of the road geometry to decide whether the vehicle will stand in its current lane or if it is doing a lane change maneuver. This is done by comparing local parabolic models of the lane's center-line and of the trajectory.

Knowing the vehicle's current state and its future lane, several possible maneuver's end points are generated by sampling in a limited area. Then, a set of trajectories are generated by joining the current vehicle's position to each end point. The trajectories are calculated so that their longitudinal and their lateral components are quintic time polynomials in the Frenet frame along the center-line of the vehicle's current lane. The MRM-based predicted trajectory is selected as the best one, with respect to a cost function that penalizes timetaking trajectories and trajectories with high pick of lateral acceleration, synonym of low comfort in the cockpit. Fig. 1 shows an example of generated trajectories and the selected one for a lane change maneuver. Even though the trajectories are defined with explicit polynomials, the conversion from the Frenet frame to the working Cartesian frame does not have an analytical solution if the shape of the road is not a straight line or a perfect circle (see appendix of [START_REF] Houenou | Vehicle trajectory prediction based on motion model and maneuver recognition[END_REF] for the conversion method). So, the trajectory is sampled and the conversion is made point by point. The MRM-based trajectory prediction (denoted T man in the following) is thus T where (x, y)

is the position in the considered Cartesian frame, θ the yaw angle, v and a the longitudinal velocity and acceleration and ω the yaw rate. Defining X 0 as the current estimated state, Eq.1 gives the closed-form of the predicted trajectory (denoted T mdl in the following). Index 0 expresses the current estimated values.

T mdl :

                     x (t) = a0 ω 2 0 cos (θ (t)) + v(t) ω0 sin (θ (t)) + cx y (t) = a0 ω 2 0 sin (θ (t)) -v(t) ω0 cos (θ (t)) + cy θ (t) = ω 0 • t + θ 0 v (t) = a 0 • t + v 0 a (t) = a 0 ω (t) = ω 0 (1)
where cx and cy are constant values obtained from the current state's values.

cx = x 0 -v0 ω0 sin (θ 0 ) -a0 ω 2 0 cos (θ 0 ) cy = y 0 + v0 ω0 cos (θ 0 ) -a0 ω 2 0 sin (θ 0 )
The final predicted trajectory, denoted T f in is a linear combination of T man and T mdl :

T f in (t) = f (t) • T mdl (t) + (1 -f (t)) • T man (t) (2)
where f (t) is a predefined function as the one depicted in Fig. 2. Fig. 3 shows an example of a trajectory prediction. One can notice that since we do not have an explicit expression of T man , T f in is also a set of predicted sampled poses. The problem is now to estimate the covariance matrix of the prediction error of each future pose.

III. PROPAGATION OF UNCERTAINTY

Let's define ξ = [x, y, θ] T as the pose of a considered vehicle at a given time. (x, y) is the position in the working Cartesian frame and θ is the yaw angle. The estimated current pose ξ 0 and its error covariance matrix Σ 0 are known from the tracking system. Let

ξ N = {ξ 1 , ξ 2 , • • • , ξ k , • • • ξ N }
be the set of predicted poses with a sampling period denoted T e . If the predicted trajectory was known in a closed-form as ξ k+1 = g (ξ k ) + α, with g (ξ) a known function and α the model noise, the covariance matrix of ξ k+1 could be calculated with a first order expansion:

Σ k+1 = J • Σ k • J T + Q (3)
where Σ k and Q are respectively the covariance matrices of the error on ξ k and of α. J = ∂g ∂ξ (ξ k ) is the Jacobian of function g in point ξ k . In our case, such a g function is not known. So, we first consider a new state vector X defined as

X = ξ ξ =         x y θ v x v y ω        
v x and v y are respectively the velocities along x-axis and yaxis, ω is the yaw rate. The error covariance matrix of X is denoted P. The initial (current) values X 0 and P 0 are known from the tracking system. Then, we assume a constant yaw rate and velocity movement between two consecutive samples X k and X k+1 , so that v x,k , for instance, is given by:

v x,k = x k+1 -x k T e
If one integrates along the trajectory with x k+1 = x k + T e .v x,k and develops this expression, one can check that we have x k+1 = x k+1 . Therefore, starting from X 0 , and having the same reasoning for all the state parameters, we get exactly the same points of the predicted trajectory. The linear state evolution is:

X k+1 = A.X k + α (4)
where

A =         1 0 0 T e 0 0 0 1 0 0 T e 0 0 0 1 0 0 T e 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 1         (5)
and the zero-mean white model noise α is added to handle the fact that this model is an approximation of the real movement (Fig. 4). In practice, the covariance matrix Q of α is chosen constant with its first three diagonal elements null because the noise does not affect the next pose (the first three state elements) which is already known from the trajectory prediction. The error covariance matrix of X k+1 is given by:

P k+1 = A • P k • A T + Q (6)
The covariance matrix Σ k+1 of ξ k+1 is finally extracted as the upper left block of P k+1 . Fig. 4 shows a simplified example (with only position and velocity) of the propagation of the initial uncertainty along a predicted trajectory. In this example:

P 0 =     1 0.81 0 0 0.81 1 0 0 0 0 1 0.64 0 0.64 0 1     , Q =     0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 1    
The model noise here, is chosen pretty large in order to easily check the increase of the uncertainty, depicted by ellipses in the figure.

IV. ESTIMATION OF THE COLLISION RISK

The risk of collision is computed separately for each obstacle in a limited area around the ego vehicle. In [START_REF] Laugier | Probabilistic analysis of dynamic scenes and collision risks assessment to improve driving safety[END_REF], the predicted trajectory of each considered object (ego vehicle or target object) is used as mean value of a Gaussian Process [START_REF] Rasmussen | Gaussian processes for machine learning[END_REF] and several random trajectory samples are generated for that object. Then, the fraction of the ego vehicle's trajectories that collide with one or more trajectories of a considered object gives the risk of collision with that object. This approach gives an idea of the collision risk on the whole predicted trajectory but does not provide an estimation of theTTC nor the point of the collision which are yet prime indicators, in the sense of collision avoidance.

In [START_REF] Berthelot | Handling uncertainties in criticality assessment[END_REF], the predicted trajectories are used to define a time function representing the distance between the ego vehicle and the considered obstacle. The smallest positive zero (if any) of this function stands for the TTC. Then, the predicted poses of both objects at that time and their variances are used to generate several Gaussian random poses. Finally, the fraction of poses that depict a collision gives the probability collision at that time. Here, the variances of the future poses are considered equal to their initial values (current time). This approach uses a deterministic method to determine the TTC, which might be a fail since the predicted trajectories are not certain. The actual TTC may potentially be different. Moreover, the variances actually increase in time.

In [START_REF] Lambert | A fast monte carlo algorithm for collision probability estimation[END_REF], the authors proposed a method to compute the probability of collision of two objects at a given moment, knowing an estimate of their respective poses at that time and the corresponding error covariance matrices. The poses are considered as Gaussian variables and the probability that the objects (named n and m) share the same space is:

P coll = ˆD p n (ξ) p m (ξ) dx n dy n dθ n dx m dy m dθ m ( 7
)
where the probability density functions p n and p m are Gaussian:

p i (ξ) = 1 √ 2π 3 |Σ i | e -1 2 ((ξ-ξi) T Σ -1 i (ξ-ξi)) |i = {n, m} and D = (xn, yn, xm, ym) ∈ R 4 , (θn, θm) ∈ [0, 2π[ 2 |Sn ∩ Sm = ∅
S n and S m being the areas respectively occupied by objects n and m. If S n = S m = ∅, Eq. 7 has an analytical solution, otherwise it does not exist. Since we consider objects with non null surface, the probability is computed numerically via a Monte-Carlo (MC) simulation, consisting in generating random poses for each object and computing the ratio of pairs of poses corresponding to a collision. In a similar way to the previously cited paper, our approach is also based on a Monte Carlo simulation but it overcomes the weaknesses mentioned before. Its main steps are the following:

1) For each considered object (including ego vehicle), propagate the initial uncertainty along its predicted trajectory (see section III) 2) For each future time sample, compute the probability of collision of the ego vehicle with each target object.

The probability is approximated by generating random poses for the ego vehicle and for the object, based on their predicted poses and associated covariance at this time. Then, the ratio of cases corresponding to a collision is the probability of collision. 3) For each target object, output the set of probability values (corresponding to the different future time samples).

Instead of estimating the actual TTC for each target object, each future time sample is considered as a potential TTC and the probability of collision is computed. This approach allows to have a larger view of the evolution of the scene and to set up different strategies, according to the probability of collision for specific values of TTC.

V. SIMULATION AND RESULTS

The risk assessment approach has been tested in simulation. Two different use cases were evaluated, with the goal of checking the evolution of the probability of collision of the ego vehicle with each encountered object. Both simulation last 12s. There are 3 vehicles, the ego vehicle named V 0 and two target vehicles named V 1 and V 2 . The curves representing the evolution of the risk of collision for the next 4s are plotted all along each simulation, as shown in Fig. 5 and Fig. 6. The upper sub-graphs shows a bird view of the scenario at the mentioned time instant. We skip the letter "V " in the name of the vehicles for a better readability of the figures. So, "0" stands for V 0 , "1" for V 1 and so on.

The middle and the bottom sub-graphs depict the risk of collision of V 0 respectively against V 1 and V 2 . Both scenarios happen on a straight road with 2 lanes, each 3.5m wide.

The simulator provides position, velocity and orientation data for each vehicle, at a rate of 10Hz. At each sample time, a 4s trajectory prediction is made for each object. Then, the predicted trajectories are sampled, also on a 10Hz rate and the probability of collision of V 0 against each vehicle is computed for each corresponding pair of future sample poses. For the MC algorithm, only 100 drawings are made per predicted pose. This value may seem low but gives realistic results as discussed in the following subsections. Moreover, it is a good compromise to afford a real-time execution. One may try a higher number of drawings to have a better approximation of the integral in Eq. 7, if the required computational resources are available.

A. Scenario 1: Collision during an overtaking

In the first scenario, the velocities of V 0 , V 1 and V 2 are respectively 80km/h, 50km/h and 40km/h. At the beginning, all the vehicles are in the bottom lane. V 0 follows V 1 which follows V 2 . Then, V 0 starts an overtaking on V 1 by making a lane change to the upper lane. While it gets close to V 1 , the latter suddenly starts a lane change in order to overtake V 2 . A collision happens in the upper lane, between V 0 and V 1 at time 7.4s. Fig. 5 shows the risk of collision at different moments of the simulation. As it can be noticed from Fig. 5, the pick of probability of collision for V 1 is always around time 7.4s (when reading the figure, please add the prediction time to the current simulation time) and the value of the pick increases as the car gets closer to the moment of the collision. This is the desired behavior because it allows a CAS to react proportionally to the time before the collision. This also means that the trajectory prediction method and the collision risk estimation are consistent.

By checking the values of the curve for V 1 , we have about 50% and 65% respectively 2s and 1.5s before the collision, which are medium levels of alert and leaves enough time to make a maneuver in order to avoid the collision. 1s before the collision, the probability reaches 80% and reaches 100% around 0.8s meaning that the collision is very soon and that there are very few chances to avoid it if no action is taken immediately.

Regarding V 2 , the risk of collision is not null despite the fact that it is not in the same lane as V 0 . Indeed, when they are predicted to be almost side by side or at least pretty close to each other, the variances of the predicted poses leave a few chances to have a collision. Nevertheless, the maximal pick is 40% and the mean value along the whole simulation is less than 20% which means that the risk of collision between V 0 and V 2 remains, as expected, low.

B. Scenario 2: Front to front collision narrowly avoided

In the second scenario, V 0 follows V 1 in the bottom lane while V 2 is driving in the upper lane but in the opposite direction. They respectively drive at 50km/h, 40km/h and 40km/h. Then, V 0 starts a lane change in order to overtake V 1 .

A collision has almost happened between V 0 and V 2 shortly later but has narrowly been avoided because V 0 came back to the bottom lane just in time. The goal here, is to check the behavior of the system when the driving behavior changes. Fig. 6 shows the risk of collision at different moments of the simulation.

In this scenario, there was no collision and, as expected, the risk of collision with V 2 never reached 100%. It even remains quite low despite the fact that both vehicle happened to be driving in opposite directions in the same lane. This is so because at the beginning, even though they were in the same lane, they are pretty far from each other (more than 150m in the first sub-figure of Figure 6).

Later, the lane change made by V 0 has been detected early making their predicted trajectories not meet. However, as the vehicles get closer, the risk of collision increases. The pick values reached in the two bottom sub-figures of Figure 6 show that they were close to a collision. Please notice that the trajectory prediction is crucial to obtain a realistic result. For V 1 the analysis is similar to the one for V 2 in the first scenario.

Figure 1 .

 1 Figure 1. MRM-based trajectory prediction

Figure 3 .

 3 Figure 3. Example of a trajectory prediction at a given time during a lane change. T mdl gets wrong as time increases.T f in is closer to the actual trajectory than each of Tman and T mdl .

Figure 4 .

 4 Figure 4. Illustration of the propagation of an initial uncertainty along a predicted trajectory.

Figure 5 .

 5 Figure 5. Probability of collision during simulation N°1

VI. CONCLUSION

This paper presented an approach to assess the risk of collision of the ego vehicle in a dynamic scene. This, first includes object trajectory predictions. So, we reminded a trajectory prediction method which was used later in the evaluation of the approach. In order to calculate the risk of a predicted collision, it is crucial to know the uncertainty of the predicted poses. A method was presented to propagate the estimated uncertainty of the current pose of a vehicle along its predicted trajectory. Then, a method to approximate the risk of collision between two objects at a given time was shown.

It is based on a Monte Carlo strategy. The global idea of the risk assessment is to have a pseudo-continuous evolution of the probability of collision on the whole prediction horizon so that it possible for a Collision Avoidance System to set up different action strategies depending on the value of the probability at specific future time instants. The overall approach was tested in two different scenarios and the results show that the output is very realistic which should allow to reduce false alarms in a CAS when the driving behavior changes before a coming collision.


Louisa Chouguiat

Raouf Boulcina

Bertrand Carboni

Albert Demonceau

Abdelmadjid Debache
email: a_debache@yahoo.fr.

A new and efficient one-pot synthesis of 2-hydroxy-1,4dihydrobenzoxazines via a three-component Petasis reaction

Keywords: 2-Hydroxy-2H-1, 4-benzoxazine, Petasis reaction, Secondary amines, Multicomponent reaction

à 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.

A new and efficient one-pot synthesis of 2-hydroxy-1,4-dihydrobenzoxazines via a three-component Petasis reaction

Graphical Abstract

A new and efficient one-pot synthesis of 2hydroxy-1,4-dihydrobenzoxazines via a three-component Petasis reaction Louisa Chouguiat, Raouf Boulcina, Bertrand Carboni, Albert Demonceau, Abdelmadjid Debache Leave this area blank for abstract info.

NH 2 OH R 1 R 4 R 3 R 2 O H + N R 5 O R 1 OH R 2 R 3 R 4 two steps H O H O R 5 B OH

OH major

Multicomponent reactions are highly efficient and atomeconomic transformations in synthetic organic chemistry. 1 They can be used for constructing an array of compound libraries in the medicinal chemistry field. The Petasis reaction (the application of a boronic acid nucleophile in the Mannich reaction) has received increasing attention due to its utility as a powerful and convenient method for the synthesis of functionalized amine derivatives, 2 such as α-amino carboxylic acids and derivatives, 3 β-amino alcohols, 4 allyl amines, 5 and various heterocyclic compounds, 6 while new applications continue to be developed. 7 This one-step, three-component reaction involves coupling of an amine, an organoboronic acid, and a carbonyl compound functionalized at the α-position to give the corresponding amine derivative in a short and experimentally simple process. Along with the accessibility of reagents and mild reaction conditions, this approach is an attractive alternative to other methodologies. A large variety of alkenyl, aryl, and heteroaryl boronic acids, various primary and secondary amine derivatives (e.g., diamines, N-hydroxylamines, N-sulfinyl amines, hydrazines) and functionalized carbonyl compounds (e.g., α-keto acids and α-hydroxy aldehydes) have been shown to participate in this reaction with success. [START_REF] Candeias | [END_REF] The reaction can be carried out in many different solvents, usually CH2Cl2, toluene, alcohols as methanol or hexafluorooisopropanol and also in water. 8 Additionally, the 1,4-benzoxazine 7 moiety is an integral part of several naturally occurring substances. For example, various glycosides of the 2-hydroxy-2H-1,4-benzoxazines skeletons have been found to occur in graminaceous plants such as maize, wheat, rye, and rice, and have been suggested to act as plant resistance factors against microbial diseases and insects. 8 The 1,4-benzoxazine moiety was also found in various antibiotics such as C-1027. 9 Generally, 1,4-benzoxazine compounds were usually synthesized via a multistep process, such as the cyclocondensation of o-aminophenols with suitable dihalo derivatives, 10 cyclocondensation of o-aminophenols with αhalo-acyl bromides followed by carbonyl group reduction with BH3, 11 and alkylation of o-nitrophenol with a haloaldehyde followed by reductive cyclization. 12 Alternatively, these 1,4-benzoxazine moieties can be prepared via ring-opening of an epoxide with ohalosulfonamides followed by cyclization 13 or by ringopening of an epoxide with o-aminophenols followed by cyclocondensation. 14 As part of our ongoing research on multicomponent reactions, 15 we have developed, a new alternative route to the conventional multistep synthesis of 2-hydroxy-1,4benzoxazines via a Petasis multicomponent reaction (Scheme 1).

Scheme 1.

The proposed strategy involves the use of 2-(benzylamino)phenols 4 as reactant. They were initially prepared in high yields by condensation of substituted 2aminophenols 1 and benzaldehydes 2 followed by reduction in a one-pot process. These key intermediates subsequently reacted with aryl boronic acids 5 and glyoxal 6 in methanol at ambient temperature to afford the target compounds 7 in good yields.

We began our research by preparing secondary amines 4. Indeed, 2-(4-methoxybenzylamino)phenol (4a) was obtained from the condensation of 2-aminophenol (1a) (1 equiv.) with benzaldehyde (2a) (1 equiv.) in dry methanol or dichloromethane at room temperature. The reaction proceeded smoothly and provided, after 24 hours, an excellent yield of the corresponding imine (3a). The addition of an excess of NaBH3CN and stirring the mixture for 24 hours at ambient temperature resulted in the formation of the desired secondary amine 4a in very high yield. This process was generalized to synthesize various substituted amines as shown in Table 1. Most of the reactions were found to proceed cleanly in about 60-90% yield of the products. The time taken for the completion of the reactions varying between 18 and 24 hours (Scheme 2).

NH 2 OH R 1 R 4 R 3 R 2 O H + N R 5 O R 1 OH R 2 R 3 R 4 two steps H O H O R 5 B OH

OH major

Our investigations on the Petasis process identified dichloromethane at room temperature as the optimum medium for the reaction of 2-(4-methoxybenzylamino)-4methylphenol 4b with phenylboronic acid 5a and glyoxal 6. These conditions provided useful yields and reasonable reaction times. Further, we observed improved results by employing methanol as the solvent, resulting in a yield of 89% of the desired 1,4-benzoxazine 7b (Table 2,entry 2). Despite the relatively long reaction times proving to be a limitation of room temperature processes, high yields of the desired product were achieved (Scheme 3). On the basis of these results, we have explored the substrate scope and limitations. A series of diverse 1,4benzoxazines bearing different substitution patterns 7a-i were prepared in moderate to excellent yields by reactions of aryl boronic acids 5 with 2-(arylmethylamino)phenols 4a-f and glyoxal 6, and all the results are listed in Table 2. [START_REF]Typical experimental procedure for the synthesis of 2-(arylamino)phenols 4: A mixture of 2-aminophenol derivative (1) (3 mmol) and aromatic aldehyde (2) (3 mmol) in MeOH (10 ml) was stirred at room temperature for 24 h (monitored by TLC). After formation of the imine 3, NaBH3CN (9 mmol) and HCl (a few drops), were added and stirring was continued for 24 h. After the reaction was complete, the mixture was poured into ice water. The resulting precipitate was filtered, washed with H2O, and dried. The resulting amine was obtained in very high purity and in excellent yield. 2-(4-Methoxybenzylamino)phenol[END_REF] It was observed that the reactions of phenylboronic acid furnished the products in higher yields than those with 4methylphenylboronic acid.

The scope of the reaction with regard to the secondary amines was also explored. Electron-donating substituents underwent the reaction with best yields (Table 2, entries 1-3). The presence of an electron-withdrawing group caused a reduction in yield (entries 4 and 5), even more significant in the case of a nitro substituent (entries 6-9). As outlined in Table 2, all these compounds were obtained as mixtures of diastereoisomers with the ratio depending on substituents in different positions of the 1,4benzoxazine core. The assignment of a relative trans configuration for the major isomer was secured by X-ray crystallographic analysis of 7f, the hydroxyl group being in an axial position in the solid state (Figure 1). 17

H N R 4 OH R 1 R 2 R 3 MeOH, r.t. 24 h R 5 B OH OH H O H O + N R 5 O R 1 OH R 2 R 3 R 4 6 H b H a 2 3 major

  Scheme 3.

  Scheme 4.

Figure 1 .

 1 Figure 1. X-ray crystal structure (ORTEP) of compound 7f (CCDC 1012140). Crystal data:C21H18N2O4, Mr = 362.37; orthorhombic, Pna21; Hall symbol: P 2c-2n; a = 12.7332 (14) Å; b = 14.2777 (14) Å; c = 19.003 (2) Å; V = 3454.8 (6) Å 3 , Z = 8.

Table 1 .

 1 Synthesis of 2-(arylmethylamino)phenols 4. a

	Scheme 2.
	Entry	R 1	R 2	R 3	R 4	2-(Arylmethylamino)phenol	Yield b (%)	M.p. (°C)
	1	H	H	H	OCH3		83	104-106
						4a
	2	CH3	H	H	OCH3		84	136-138
						4b
	3	Cl	H	H	OCH3		60	148-150
						4c
	4	CH3	H	H	Br		90	86-90
						4d
	5	H	NO2	H	H		75	90-92
						4e
	6	CH3	H	NO2	H		95	120-122
						4f

a Reactions conditions: 1 (1 mmol), 2 (1 mmol), MeOH (5 ml), r.t., 24 h, then NaBH3CN (3 mmol) was added and stirring was continued for 24 h. b Isolated yields.

Table 2 .

 2 The Petasis reaction of secondary amines 4, aryl boronic acids 5, and glyoxal (6). a may be useful compounds for medicinal and materials chemistry. Simple substrates, good yields, and operational simplicity are significant advantages of this procedure. Further investigations regarding the exact mechanism of this Petasis condensation and the further transformation of these heterocycles are currently in progress.

	Entry	Secondary amine	R 5	1,4-Benzoxazine	M.p. (° C)	Yield b (%)	dr c
	1	4a	CH3		oil	76	95:5
				7a
	2	4b	H		oil	89	86:14
				7b
	3	4c	H		oil	70	88:12
				7c
	4	4d	H		131-132	72	90:10
				7d

Acknowledgment

The authors gratefully acknowledge le Ministère de l'Enseignement Supérieur et de la Recherché Scientifique (Algeria) for the financial support of this research program.

Supplementary data

Supplementary data associated with this article can be found, in the online version, at http://dx.doi.org/10.1016/j.tetlet.2014.07.093.