Pierre-Louis Antonsanti
email: pierrelouis.antonsanti@ge.com

Thomas Benseghir
email: thomas.benseghir@ge.com

Vincent Jugnon
email: vincent.jugnon@ge.com

Joan Glaunès

Database Annotation with few Examples: An Atlas-based Framework using Diffeomorphic Registration of 3D Trees

Automatic annotation of anatomical structures can help simplify workflow during interventions in numerous clinical applications but usually involves a large amount of annotated data. The complexity of the labeling task, together with the lack of representative data, slows down the development of robust solutions. In this paper, we propose a solution requiring very few annotated cases to label 3D pelvic arterial trees of patients with benign prostatic hyperplasia. We take advantage of Large Deformation Diffeomorphic Metric Mapping (LD-DMM) to perform registration based on meaningful deformations from which we build an atlas. Branch pairing is then computed from the atlas to new cases using optimal transport to ensure one-to-one correspondence during the labeling process. To tackle topological variations in the tree, which usually degrades the performance of atlas-based techniques, we propose a simple bottom-up label assignment adapted to the pelvic anatomy. The proposed method achieves 97.6% labeling precision with only 5 cases for training, while in comparison learningbased methods only reach 82.2% on such small training sets.

Introduction

The automatic annotation of tree-like structures has many clinical applications, from workflow simplification in cardiac disease diagnosis ( [START_REF] Akinyemi | Automatic labelling of coronary arteries[END_REF], [5], [START_REF] Ezquerra | Model-guided labeling of coronary structure[END_REF], [START_REF] Wu | Automated anatomical labeling of coronary arteries via bidirectional tree lstms[END_REF]) to intervention planning in arterio-venous malformations ( [START_REF] Bogunović | Anatomical labeling of the circle of willis using maximum a posteriori probability estimation[END_REF], [START_REF] Ghanavati | Automatic anatomical labeling of the complete cerebral vasculature in mouse models[END_REF], [START_REF] Robben | Simultaneous segmentation and anatomical labeling of the cerebral vasculature[END_REF], [START_REF] Wang | Automatic labeling of vascular structures with topological constraints via hmm[END_REF]) and lesion detection in pneumology ( [START_REF] Feragen | Geodesic atlas-based labeling of anatomical trees: Application and evaluation on airways extracted from ct[END_REF], [START_REF] Van Ginneken | Robust segmentation and anatomical labeling of the airway tree from thoracic ct scans[END_REF], [START_REF] Lo | Vascular image registration techniques: A living review[END_REF]). This task is particularly important in the context of interventional radiology, where minimally invasive procedures are performed by navigating small tools inside the patient's arteries under X-ray guidance. For example, benign prostatic hyperplasia symptoms can be reduced by embolizing arteries feeding the prostate to reduce its size [START_REF] Ray | Efficacy and safety of prostate artery embolization for benign prostatic hyperplasia: an observational study and propensity-matched comparison with transurethral resection of the prostate (the uk-rope study)[END_REF].

Identifying the correct arteries to treat during the procedure -along with their neighbours -is crucial for the safety of the patient and the effectiveness of the treatment. However, this task is very challenging because of the arterial tree complexity and its topological changes induced by frequent anatomical variations. Having a large representative annotated database is also a challenge, especially in medical imaging where the sensitivity of the data and the difficulty to annotate slow down the development of learning based techniques [START_REF] Lee | Medical big data: promise and challenges[END_REF]. In this context, solutions to the tree labeling problem should ideally work from a limited number of annotated samples.

Learning-based labeling With the prevalence of machine learning, many articles address the problem of automatic anatomical tree labeling by extracting features from the tree to feed a learning algorithm that predicts labels probabilities. Authors of [START_REF] Van Ginneken | Robust segmentation and anatomical labeling of the airway tree from thoracic ct scans[END_REF] and [START_REF] Akinyemi | Automatic labelling of coronary arteries[END_REF] first proposed Gaussian mixtures models to learn from branch features (geometrical and topological) and predict the label probabilities. In [START_REF] Akinyemi | Automatic labelling of coronary arteries[END_REF], labels are assigned following clinical a priori defining a set of rules. In [START_REF] Lo | Vascular image registration techniques: A living review[END_REF], k-Nearest-Neighbours predicts label probabilities that are used in a bottom-up assignment procedure, searching through all the existing branch relationships in the training data. Since numerous features leads to high dimensional space, such techniques often show poor generalization capacity.

To reduce dimensionality, boosting algorithms and tree classifiers are used in [START_REF] Hoang | A study on automated anatomical labeling to arteries concerning with colon from 3d abdominal ct images 7962[END_REF], [START_REF] Matsuzaki | Automated anatomical labeling of abdominal arteries and hepatic portal system extracted from abdominal ct volumes[END_REF] and [START_REF] Robben | Simultaneous segmentation and anatomical labeling of the cerebral vasculature[END_REF] to select the most discriminant features. Assignment procedure is enriched by topological rules to improve coherence along the tree, instead of considering branches independently. Label assignment was further refined in [START_REF] Ghanavati | Automatic anatomical labeling of the complete cerebral vasculature in mouse models[END_REF] and [START_REF] Wang | Automatic labeling of vascular structures with topological constraints via hmm[END_REF] by adding a Markovian property -each branch label depending on its direct neighbours -with Markov chains parameters learned during training. Such assignment strongly depends on the database size and the task complexity in term of anatomical variability.

While previous articles had access to limited size databases (around 50 cases), authors of [START_REF] Wu | Automated anatomical labeling of coronary arteries via bidirectional tree lstms[END_REF] trained a recurrent neural network preserving the topology of the coronary tree on a database composed of 436 annotated trees. The labels predictions rely on a multi-layer perceptron and a bidirectional tree-structural long short-term memory network. This interesting approach using deep learning for vessel classification is less common because it requires a lot more training data to be able to capture the anatomical variability.

Atlas-based labeling Contrary to learning-based techniques, atlas-based ones can offer robust annotation even with few annotated cases ( [START_REF] Wu | Automated anatomical labeling of coronary arteries via bidirectional tree lstms[END_REF]). An atlas is defined as a reference model that can be built from prior knowledge ([3], [5]) or from an available annotated database ( [START_REF] Bülow | Point based methods for automatic bronchial tree matching and labeling 6143[END_REF], [START_REF] Ezquerra | Model-guided labeling of coronary structure[END_REF], [START_REF] Feragen | Geodesic atlas-based labeling of anatomical trees: Application and evaluation on airways extracted from ct[END_REF]). Most of the atlas-based methods follow a four steps framework: the choice of the atlas, the registration onto the target, the estimation of the labels probabilities and finally the assignment. In [START_REF] Ezquerra | Model-guided labeling of coronary structure[END_REF], authors focused on the annotation procedure, relying on previously established atlas and a manual registration. Labeling is done through a branch-and-bound algorithm extending the best partial labeling with respect to a function designed to compare observations to the atlas.

Later in [START_REF] Bülow | Point based methods for automatic bronchial tree matching and labeling 6143[END_REF], label probabilities are computed at the level of branches by a voting procedure, each point along a branch voting for a label. This solution can adapt to missing branches but still does not guarantee anatomical consistency of the labeling along the tree and is not robust to topological variation with respect to the atlas. To take this variability into account [START_REF] Bogunović | Anatomical labeling of the circle of willis using maximum a posteriori probability estimation[END_REF] proposes to create one atlas per known topology in the anatomy of interest, an interesting approach if the number of anatomical variants is limited. Similarly [5], [START_REF] Feragen | Geodesic atlas-based labeling of anatomical trees: Application and evaluation on airways extracted from ct[END_REF] use multi-atlas approaches, taking advantage of a distance to the atlas that quantifies topological differences. In [5], the reference case is selected as the best example in a training set following a leave-one-in cross-validation design.

While a lot of advanced methods have been proposed for medical structure registration [START_REF] Sotiras | Deformable medical image registration: A survey[END_REF], few efforts have been made to use advanced deformation models in the context of vascular tree labeling [15]. In our work, we propose an atlas-based algorithm Fig. 1: The proposed atlas-based vascular tree annotation pipeline. illustrated in Fig. 1 relying on state of the art LDDMM to build the atlas and estimate realistic deformations. It is combined with an Optimal Transport based matching, to compute a relevant assignment between atlas and target branches. In the end, to handle the topological changes between the atlas and the target, a bottom-up label assignment is performed to achieve optimal results for our pelvic vasculature labeling problem.

Method

Let a labeled vascular tree T = ({(B α , l α )} α , M ) be a set of branches B α labeled l α ∈ N with connections to other branches stored in an adjacency matrix M . Two labeled branches (B α , l α ) and

(B β , l β ) are connected if M α,β = 1. Each branch B α is a polygonal curve composed of ordered points B α [k] ∈ R 3 representing a vessel centerline.
We also denote {q i } i the unordered set of all points B α [k] of the tree for all indices α and k.

The Large Deformation Diffeomorphic Metric Mapping framework We compute the deformation ϕ of a source shape S (in our case the atlas, a centerlines tree) onto a target shape T using LDDMM. This state of the art framework, detailed in [START_REF] Younes | Shapes and Diffeomorphisms[END_REF] (chapters 8 to 11), allows to analyze differences between shapes via the estimation of invertible deformations of the ambient space that act upon them. In practice, the diffeomorphism ϕ is estimated by minimizing a cost function J (ϕ) = E(ϕ) +A(ϕ(S), T ) where E is the deformation cost, and A is a data attachment term that penalizes mismatch between the deformed source ϕ(S) and the target T . We chose the Normal Cycles model proposed in [START_REF] Roussillon | Kernel metrics on normal cycles and application to curve matching[END_REF] as data attachment since it has shown good results at areas of high curvature and singular points as bifurcations or curves endpoints.

The deformation map ϕ : R 3 → R 3 is defined as the solution at time t = 1 of a flow equation ∂ t φ(t, x) = v(t, φ(t, x)) with initial condition φ(0, x) = x, where v(t, •) are time varying vector fields assumed to belong to a Hilbert space V , ensuring regularity properties. As shown in [START_REF] Miller | Geodesic shooting for computational anatomy[END_REF], in a discrete setting, denoting x i , i ∈ [1, ..., n S ] the discretization points of the source shape, one may derive optimality equations that must be satisfied by the trajectories q i (t) = φ(t, x i ) when considering deformations that minimize the cost function. These equations take the following form:

             ṗi (t) = - 1 2 ∇ qi(t)   n S j=1 n S l=1 p j (t) , K V (q j (t), q l (t))p l (t)   qi (t) = n S j=1 K V (q i (t), q j (t))p j (t) (1)
where p i (t) are auxiliary dual variables called momenta. They correspond to geodesic equations with respect to a specific Riemannian metric, written in Hamiltonian form.

The deformations are then fully parametrized by the set of initial momenta

(p i (0)) i ∈ R 3n S .
To model deformations, we define the kernel to be a sum of Gaussian kernels

K V (x, y) = s exp -x -y 2 / (σ 0 /s) 2
, where s ∈ [START_REF] Akinyemi | Automatic labelling of coronary arteries[END_REF][START_REF] Bülow | Point based methods for automatic bronchial tree matching and labeling 6143[END_REF][START_REF] Feragen | Geodesic atlas-based labeling of anatomical trees: Application and evaluation on airways extracted from ct[END_REF][START_REF] Matsuzaki | Automated anatomical labeling of abdominal arteries and hepatic portal system extracted from abdominal ct volumes[END_REF]] and σ 0 is half the size of our vascular trees bounding box. This multi-scale approach was introduced by [START_REF] Risser | Simultaneous multi-scale registration using large deformation diffeomorphic metric mapping[END_REF]. Normal Cycles also require to choose specific kernels; we take a sum of constant and linear kernels for the spherical part and a Gaussian kernel for the spatial part, at two scales σ 0 and σ 0 /4, using the output (p i (0)) i of the optimization process at scale σ 0 as initialization of the optimization at scale σ 0 /4.

Building the atlas In the context of vascular tree labeling, the tree topology can be highly variable (in our database we have one topology every two cases). Trying to build atlases that take into account these variations is unrealistic. Yet we still want the automatic annotation to be robust to the choice of the reference case.

The LDDMM framework allows to derive methods for computing such atlases. More precisely, the optimal deformations generated by LDDMM are fully parameterized by the initial momenta, and their representation in the Euclidean space R 3n S allows to perform classical linear statistics. Following [START_REF] Vaillant | Statistics on diffeomorphisms via tangent space representations[END_REF], in order to build the atlas we select one available annotated case and compute its deformations onto a set of N targets (the selected case included). These registrations provide a collection of initial momenta

{ p k i (0) i , k ∈ [1, ..., N ]}.
The reference case is selected using a leave-one-in method as in [5] with respect to the labeling procedure. The atlas is then obtained by shooting via geodesic equations to deform the selected case along the average of the initial momenta p i (0) = (1/N ) k p k i (0). This procedure does not require the targets to be annotated, and seems suited to build an atlas of the whole database.

In addition it can be iterated by replacing the selected case with the deformed one. We refer to LDDMM-k for the k-th iteration of the atlas construction. Consequently LDDMM-0 refers to the case of using the reference case directly as an atlas. As illustrated in Fig. 2 the atlas converges along the iterations to an average position representative of the set of targets. Through this spatial normalization, we limit the sensitivity to the choice of the initial case as atlas. To illustrate the impact of LDDMM-k on the labeling, we chose a simple assignment procedure described in [START_REF] Bülow | Point based methods for automatic bronchial tree matching and labeling 6143[END_REF]: each point of the target tree votes for the label of its closest point in the labeled atlas. Then, we compute the vote per branch B T β in the target:

π(B T β , l) = Card({q T i ∈ B T β , l(q T i ) = l}) / Card(B T β
). This label probability estimation does not guarantee anatomical consistency and deeply relies on the registration. Additionally, each point vote is independent from the others, allowing to characterize the quality of the registration from a labeling point of view.

Optimal Transport for a better assignment During the LDDMM atlas construction and registration, each tree is seen as one shape. Consequently, there is no assumption over branch matching and topological changes. In order to provide a relevant label assignment that takes the mutual information into account we propose to use Optimal Transport. It is convenient to compute the optimal one-to-one assignment between branches of the deformed source and the target with respect to a given distance.

Based on the work of [START_REF] Feydy | Fast and Scalable Optimal Transport for Brain Tractograms[END_REF], each branch is re-sampled with 20 points and the distance matrix D between each branch B ϕ(S) α and B T β is given by:

D α,β = B ϕ(S) α -B T β R 3d .
We tried different numbers of points d per branch ranging from 20 to 500 with no significant impact over the matching results. Considering that our problem is of limited size (17 branches per tree), a simple Kuhn-Munkres algorithm (also called Hungarian algorithm) was used to compute the assignment solution. It consists in finding minimum weight matching in bipartite graphs by minimizing the function α,β D α,β .X α,β with (X α,β ) ∈ {0, 1} 17•17 the output boolean matrix with 1 if the branches B ϕ(S) α is assigned to B T β , 0 otherwise. To be consistent with the alternative simple label probability estimation of our LDDMM-k pipeline, we similarly define here π(B T β , l S α ) := X α,β , although these "probabilities" are always 0 or 1 in this case. This assignment process is complementary to the LDDMM-k process since it focuses on assignment between branches while LDDMM-k focuses on the atlas construction and the registration. We will call this pipeline LDDMM-k+OT. We will also experiment the Optimal Transport assignment without any registration (i.e. taking ϕ = id), which we denote OT.

Label assignment post-processing A first label assignment procedure directly takes the highest label probability for each branch: l(B T α ) = arg max l (π(B T α , l)), π being the output of the OT procedure or the voting probabilities estimation. This direct assignment is not based on a priori knowledge and directly reflects the performance of the prediction methods.

In practice, the vessels are labeled by the expert accordingly to the anatomy or area they irrigate [START_REF] Assis | Pelvic arterial anatomy relevant to prostatic artery embolisation and proposal for angiographic classification[END_REF]. In the application to pelvic vascular tree, when two branches of different labels share a parent, this parent is called "Common Artery". This is the only clinical a priori we introduce in the method. To limit the effect of the topological variations between the atlas and the target, we propose a bottom-up assignment procedure: first for all B T α leaf of T , l(B T α ) = arg max l (π(B T α , l)) then recursively every parent branch B T α is assigned a label with the rule:

l(B T α ) = l if l(B T β ) = l for every branch B T β child of B T α 0 (i.e."Common Artery"), otherwise.
This recursive assignment procedure, although specific to this anatomy, is quite adaptable. In fact, in most of the structured tree-shaped anatomies (coronary [START_REF] Akinyemi | Automatic labelling of coronary arteries[END_REF], [5], airway tree [START_REF] Lo | Vascular image registration techniques: A living review[END_REF], pelvic [START_REF] Assis | Pelvic arterial anatomy relevant to prostatic artery embolisation and proposal for angiographic classification[END_REF]) the branches names are also conditioned by the leaves labels. When two arteries of different labels share a common parent, this parent is either unnamed (as in our application), or named by a convention provided by the experts. The latter situation corresponds to additional conditions (as in [START_REF] Lo | Vascular image registration techniques: A living review[END_REF]) during the assignment.

Results

We conducted experiments on a dataset of 50 pelvic vascular trees corresponding to 43 different patients, some trees being the left and right vasculature of a single patient. The centerlines composing the vascular tree are constructed from 3D volumes (injected Cone Beam Computed Tomography). While the entire vascular tree is composed of up to 300 different branches, we manually extracted a simplified tree composed of the main arteries documented in the literature [START_REF] Assis | Pelvic arterial anatomy relevant to prostatic artery embolisation and proposal for angiographic classification[END_REF]. This allowed us to reduce the problem to the annotation of a 17-branches binary tree that corresponds to the typical size of trees found in the literature [5], [START_REF] Wang | Automatic labeling of vascular structures with topological constraints via hmm[END_REF]. The selected arteries of interest are: the prostatic, the superior vesicle, the obturator (2), the pudendal, the inferior gluteals (2) and the superior gluteals [START_REF] Assis | Pelvic arterial anatomy relevant to prostatic artery embolisation and proposal for angiographic classification[END_REF]. Unlabeled branches in the tree correspond to proximal common portions of the arterial tree that we label "common" arteries for a total of 7 labels. This simplified representation still captures the anatomical variability described in [START_REF] Assis | Pelvic arterial anatomy relevant to prostatic artery embolisation and proposal for angiographic classification[END_REF], as we found 28 different tree labels arrangements among the 50 cases with high variability of branch shapes and positions as is illustrated in Fig. 2. Registrations are computed with the library KeOps [START_REF] Charlier | Kernel operations on the gpu, with autodiff, without memory overflows[END_REF] allowing fast GPU computing and automatic differentiation through PyTorch ( [START_REF] Paszke | Pytorch: An imperative style, highperformance deep learning library[END_REF]) backends. The optimization of the functional J is performed using Limited-memory Broyden-Fletcher-Goldfarb-Shanno (L-BFGS) algorithm. In Fig. 2 (c) we plot the precision of the LDDMM-k at each iteration k.

A first experiment was carried out to illustrate the contributions of the atlas construction on the labeling of the database. We computed LDDMM-0 (pure registration) and LDDMM-1 (atlas construction) by registering each of the 50 available cases onto the others. It is illustrated in Fig. 2 (a,b). The average precision of one reference case used in LDDMM-0 to annotate the 49 other trees is 93.3%(±3.5) when associated to bottom-up assignment and 84.2%(±4.4) using direct assignment. This drop of performance illustrates the sensitivity of atlas-based methods to the choice of the atlas in the first place. The bottom-up assignment post processing allows to overcome this sensitivity: we use it in the rest of the experiments. We then select one of the worst cases in the database regarding LDDMM-0 labeling performance and iteratively build the new atlas following the LDDMM-k procedure. We can see that performance improves with iterations, which indicates that the atlas gradually captures the database variability: it allows a better registration hence a better label assignment. This single-case solution allows to annotate the 49 cases of the database with a precision reaching 98.9%(±0.33) while being one of the worst with LDDMM-0. It must be pointed out that the atlas construction did not rely on any other annotated case than the one initially selected. We also computed LDDMM-1 using each tree as reference case and obtained an average score of 96.8%(±0.34), improving the average precision by 3.5%. In addition the standard deviation for LDDMM-1 confirms that iterating to build an atlas makes the LDDMM-k method robust to the initial choice of the atlas. In Fig. 2 (c) we observe that one iteration of the atlas construction is enough to greatly improve the labeling of the entire database, then the performance slowly increases until iteration 4. Therefore, and to limit the computational cost, we chose to use LDDMM-1 for the rest of the experiments.

To demonstrate that atlas based techniques described in Sec. 2 perform well compared to learning based ones in case of small size database, we implemented two classification algorithms working on branch features, KNN and XGboost, inspired from the work of [START_REF] Lo | Vascular image registration techniques: A living review[END_REF] and [START_REF] Wang | Automatic labeling of vascular structures with topological constraints via hmm[END_REF]. To be close to the work of [START_REF] Lo | Vascular image registration techniques: A living review[END_REF] and [START_REF] Matsuzaki | Automated anatomical labeling of abdominal arteries and hepatic portal system extracted from abdominal ct volumes[END_REF], each branch is represented by a vector composed of 28 branch features such as direction, length or geometrical characteristics of the centerline curves, and 13 tree features involving branch's relationships to the root and to its children to introduce topological information. These features are listed in the complementary material. Tuning the XGBoost parameters had no impact so we kept parameters provided in [START_REF] Wang | Automatic labeling of vascular structures with topological constraints via hmm[END_REF]. We performed 8-neigbours KNN in the space of the 10 most informative features selected by XGboost.

To compare all methods, precision has been evaluated using a cross-validation over the 50 cases with a training set and test set of varying size: from 2% (1 case) to 80% (40 cases) of the total dataset. While the notion of cross-validation is well defined for learning-based techniques, in the case of atlas-based methods, it follows the leave-onein procedure described in Sec. 2 to select the best atlas among the training data. This case is then used to annotate the test set.

For each method we used the assignment technique giving the highest precision (direct assignment for training-based methods and bottom-up for atlas-based ones). Results are presented in Fig. 3. As expected, under 20% of training data (10 cases) precision of learning-based methods drastically drops. KNN is outperformed by other approaches, and XGBoost seems to asymptotically reach the atlas-based performances. On the other hand the atlas-based methods with bottom-up assignment perform with very little influence of the size of the training set.

Despite topological variations, the LDDMM-1 approach generates meaningful registrations showing good results when coupled with bottom-up assignment. OT also gives relevant branch matching that provides the same level of performance with the bottom-up assignment. Consequently LDDMM-1+OT with the bottom-up assignment have the best results, particularly in the case of small training sets. We illustrate in Fig. 3 (b) the performances of atlas-based methods for only 5 cases in the training set (Confusion matrices for this setting are provided in the complementary material). The results of LDDMM-1+OT are significantly better than each method taken independently with an average 97.6%(±0.97) precision.

Conclusion and perspectives

We have proposed an atlas-based labeling method allowing to annotate a database with a 97.6% average precision using only 5 cases as training data. This level of precision isn't achieved by learning-based approaches even with 8 times more training data. Our method takes advantage of the LDDMM realistic deformations to build a meaningful atlas and register it onto the cases to annotate. Optimal Transport computes a global branch matching that, combined with a bottom-up label assignment, provides an anatomically consistent labeling. The bottom-up assignment procedure allows to tackle the anatomical variations that usually degrades atlas-based methods precision. This procedure may be specific to pelvic anatomy, however we believe it could be easily adapted to other ones ([1], [START_REF] Assis | Pelvic arterial anatomy relevant to prostatic artery embolisation and proposal for angiographic classification[END_REF], [START_REF] Lo | Vascular image registration techniques: A living review[END_REF]). In addition, the LDDMM framework allows to perform statistics over deformations that could be exploited to generate new realistic data or to detect anomalies. In further work we would like to extend our method to more complex cases such as missing labels or branches. In this perspective a promising lead would be to take topological changes into account as in [START_REF] Feragen | Geodesic atlas-based labeling of anatomical trees: Application and evaluation on airways extracted from ct[END_REF] for both the deformations and the data attachment.

Fig. 2 :

 2 Fig. 2: Building atlases via LDDMM registration. The colors represent the ground truth labels. (a) The initial trees; (b) Reference atlases at LDDMM-1 (c) Precision of one atlas over LDDMM-k iterations.

Fig. 3 :

 3 Fig. 3: (a) Comparison of the performance according to the training set size. (b) Box plot at training size 10% (5 cases).

Acknowledgement We would like to thank Arthur Rocha, engineer at Sao Paulo Hospital das Clinicas for his help and expertise on the annotation of pelvic vascular trees.


domain: phys spi

text:
Einar Knudsen
email: einar.knudsen@acoem.com

Julien Perchoux
email: julien.perchoux@laas.fr

Thierry Mazoyer
email: thierry.mazoyer@acoem.com

Francis Jayat
email: francis.jayat@laas.fr

Clément Tronche
email: clement.tronche@laas.fr

Thierry Bosch
email: thierry.bosch@laas.fr

Lower detection limit of the acousto-optic effect using Optical Feedback Interferometry 1 st

Keywords: Optical feedback, interferometry, self-mixing, acousto-optic effect

Measurement and 3D imaging of acoustic waves through the acousto-optic effect has recently been demonstrated by means of Optical Feedback Interferometry (OFI). In this paper we study experimentally the lower limits of detection of an acoustic wave using an OFI sensor. We show that the OFI sensor exhibits a linear response to acoustic power variations, and we obtain a lower limit of detection of 83 dB rms for a planar acoustic wave at 3 kHz. We also determine the equivalent displacement, that is seen by the OFI sensor at this pressure level, to be 96 pm. A deeper understanding of the limits of the technology and the quantification of the acousto-optic effect shall help improve the applications already created for the measurement of acoustic pressure waves using OFI.

I. INTRODUCTION

Optical Feedback Interferometry (OFI) is a non-contact, high-precision sensing scheme which is applied to many domains of metrology. Today we find OFI in applications such as flowmetry, absolute distance measurements, vibrometry and velocimetry [START_REF] Taimre | Laser feedback interferometry: a tutorial on the self-mixing effect for coherent sensing[END_REF]- [START_REF] Bertling | Imaging of acoustic fields using optical feedback interferometry[END_REF]. Compared to other non-contact sensing technologies, OFI is self-aligned, low cost and robust. Recent applications have applied the sensing scheme to the imaging of acoustic waves through the acousto-optic effect. Bertling et al [START_REF] Bertling | Imaging of acoustic fields using optical feedback interferometry[END_REF] and Urgiles et al [START_REF] Ortiz | Visualization of an acoustic stationary wave by optical feedback interferometry[END_REF] create visual representations of standing and propagating sound waves in air, while Dubois et al proposed an exciting application to image pressure waves in acoustic metamaterials [START_REF] Dubois | Acoustic flat lensing using an indefinite medium[END_REF].

In these applications they measure the minute changes in the optical path between the laser and the target, through the modulation of the refractive index of air caused by an acoustic pressure wave. These changes in optical path are usually far below the half wavelength of the laser and close to the limitations of the technology. One particular challenge that has not yet been treated is the quantification of the air-pressure change in an acoustic wave using this method. OFI operates by backscattering the laser light off a target, inducing coherent mixing of the returning light with the laser inner cavity field. The laser wave phase and amplitude is perturbed, and these perturbations are measured by the laser's in-package photodiode. The signal provides us with information about the target the laser is pointing on. If the target is in movement, interferometric fringes are observed. A fringe appears when the target has moved a distance of half the laser's wavelength λ [START_REF] Taimre | Laser feedback interferometry: a tutorial on the self-mixing effect for coherent sensing[END_REF], [START_REF] Veng | Fringe Disappearance in Self-Mixing Interferometry Laser Sensors: Model and Application to the Absolute Distance Measurement Scheme[END_REF], [START_REF] Lang | External optical feedback effects on semiconductor injection laser properties[END_REF].

Measurements inferior to the half wavelength of the laser has been reported. Currently the state of the art using the OFI sensing scheme is measurements of 100 pm [START_REF] Giuliani | Self-mixing laser diode vibrometer[END_REF], [START_REF] Kossl | Basilar membrane resonance in the cochlea of the mustached bat[END_REF], [START_REF] Lukashkin | A selfmixing laser-diode interferometer for measuring basilar membrane vibrations without opening the cochlea[END_REF]. Through the acousto-optic effect we measure very small changes in the optical path through the modulation of the refractive index of air in the laser beam's path. The effect has been depicted [START_REF] Bertling | Imaging of acoustic fields using optical feedback interferometry[END_REF], [START_REF] Ortiz | Visualization of an acoustic stationary wave by optical feedback interferometry[END_REF] but there are no evidence of quantification of the acoustic amplitude or the effect on the refractive index so far.

The refractive index of air influences the time it takes for light to travel from the laser to the target and back. P. E. Ciddor [START_REF] Ciddor | Refractive index of air: new equations for the visible and near infrared[END_REF] made an an exhaustive model to precisely calculate the refractive index of light under given conditions. His model operates in the visible and near infrared, even for small pressure variations.

In this paper we study the dynamic range of an OFI system measuring acoustic signals using the acousto-optic effect. Furthermore, we present a model for the conversion of measured acoustic pressure through the acousto-optic effect to an equivalent displacement in order to quantify the phenomenon. In earlier studies [START_REF] Bertling | Imaging of acoustic fields using optical feedback interferometry[END_REF], [START_REF] Ortiz | Visualization of an acoustic stationary wave by optical feedback interferometry[END_REF], pressure and frequency are controlled, but the waves undergo a free field spherical propagation. To quantify the acousto-optic effect, we have chosen to work with planar waves, where the advantages of such waves are described in sections II-A and II-B.

II. THEORY

A. OFI model and governing equations

The laser power under external feedback P F is calculated through the well established equation

P F = P 0 [1 + m cos(ω F τ ext )] , (1)
where P 0 is the free running laser power, m is a modulation index that depends on the target reinjection efficiency, and ω F is the laser frequency under feedback [START_REF] Taimre | Laser feedback interferometry: a tutorial on the self-mixing effect for coherent sensing[END_REF], [START_REF] Veng | Fringe Disappearance in Self-Mixing Interferometry Laser Sensors: Model and Application to the Absolute Distance Measurement Scheme[END_REF], [START_REF] Ortiz | Visualization of an acoustic stationary wave by optical feedback interferometry[END_REF]. An acoustic wave is a variation of air pressure with a subsequent change in refraction index n [START_REF] Ciddor | Refractive index of air: new equations for the visible and near infrared[END_REF]. Considering a change in n, the photon round-trip time in the external cavity τ ext is expressed as

τ ext = τ 0 + δτ (2)
with τ 0 = 2n 0 L ext /c where c is the speed of light, n 0 is the refractive index of air without acoustic excitation and

δτ = Lext 0 2δn(z) c dz , ( 3
)
where z is the coordinate along the laser beam. If we consider a plane pressure wave propagating in a perpendicular direction to the laser beam, the refractive index is constant along the z axis and we can express (3) as

δτ (t) = 2L ext δn(t)/c (4)
where δn(t) = δn max sin(ω a t), with ω a = 2πf a where f a is the acoustic frequency, and δn max is the amplitude of the refractive index modulation of air in a given acoustic wave.

B. Refracive index of air

The round-trip time τ ext in (1) is modulated through the change in n which is in turn modulated by the change of density of air under the effect of the acoustic pressure wave. To calculate variations of n we use Ciddor's model [START_REF] Ciddor | Refractive index of air: new equations for the visible and near infrared[END_REF]. The model includes wavelength, pressure, temperature, relative humidity and CO 2 concentration. Since the acoustooptic induced change in the optical path produces the same effect as a change in physical path, to quantify the acoustooptic effect we calculate an equivalent displacement from the change in refractive index

δd eq = L ext ∂n ∂p λ,P0,T0,RH0,RCO2,0 × δp (5)
where ∂n/∂p is the change in refractive index with regards to the change in pressure, and δp is the pressure difference in an acoustic wave, measured in Pa. λ is the laser wavelength, P 0 the initial pressure, T 0 the temperature in Celsius, RH the relative humidity and R CO2,0 the CO 2 concentration. The δn max following ( 4) is calculated for a pressure change δp using Ciddor's model, and the constants ∂n and ∂p are computed. By considering planar waves, the pressure p and hence the refractive index n is evenly distributed through the entire optical axis of the laser.

C. Experiment

A Distributed Feedback laser diode (1310P5DFB) emitting a transverse and longitudinal single mode at 1310 nm is mounted on an Arroyo Instruments laser mount. The mount includes a temperature control that is driven by a ThorLabs TED 200C. The laser beam is colimated by a ThorLabs C230TMD-C lens and aimed at a Thorlabs PF05-03-P01 mirror mounted in a Thorlabs KAD12F adjustable mount, 200 mm from the laser source. The mirror is aligned so that the beam is re-injected into the laser with the best achievable efficiency.

In order to assure planar waves, the acoustic waves are injected into a cylinder of diameter 55 mm and length 0.5 m. The acoustic wave is generated by a Visatron FRS 8 speaker, driven by a Renkforce SAP-702 amplifier that is mounted inside an acoustic impedance adapter of length 130 mm at the end of the waveguide. The interferometric system if fixed on the waveguide at 200 mm from the inlet. A G.R.A.S 46AE metrologic microphone is inserted into the waveguide over the laser beam, flush with the inner wall as to not perturb the pressure wave. The acoustic signal is acquired by an ACOEM 04dB acquisition card connected to the microphone.

The ensemble laser -waveguide -mirror, is placed on vibration absorbing foam to reduce potential vibrations propagating trough the optical table. The impedance adapter and speaker is attached to the waveguide through absorbent foam to reduce vibrations due to the mechanical movement of the speaker membrane. Two PCB 352C18 accelerometers are fixed on the distant sensor parts: one on the laser mount, the other on the mirror mount.

Once acquired, the accelerometer signals are integrated two times before being subtracted from one another to compute the relative displacement.

To determine the linearity and minimal detection level of the system, a continuous sinusoidal acoustic signal with frequency 3 kHz is sent through the waveguide and the OFI signal is recorded for signal amplitudes ranging from 70 dB rms

D. Numerical simulations

Using ( 5) we can calculate an equivalent displacement for an acoustic pressure wave under a given set of conditions. With an OREGON WMR89 weather station we measure the temperature to be 20 • C and the relative humidity to be 35%. The atmospheric pressure is 101.325 kPa, the laser's wavelength is 1310 nm and the CO 2 quantity is set to 450 ppm as proposed in [START_REF] Ciddor | Refractive index of air: new equations for the visible and near infrared[END_REF]. We calculate the refractive index n over the pressure range of a class one metrological sonometer (65 kPa to 110 kPa) with mentioned parameters as input, and we derivate to obtain ∂n for use in [START_REF] Veng | Fringe Disappearance in Self-Mixing Interferometry Laser Sensors: Model and Application to the Absolute Distance Measurement Scheme[END_REF]. Using the same input parameters we can simulate the variation of pressure in an acoustic wave and calculate the equivalent displacement for any given acoustic power.

III. DISCUSSION

In Fig. 2 we observe the system's response to acoustic pressure levels. When we decrease the acoustic power from 70 dB rms to 112 dB rms the response is fairly linear before flattening out at the bottom. The system's response to acoustic the pressures 87 dB rms -92 dB rms remains at a similar level. This is possibly caused by a change in ω F (1), due to a perturbation in the laser current, change in laser temperature, or to vibrational perturbations. The lowest detectable acoustic amplitude is 83 dB rms which is recorded with an amplitude of 2.4 mV by the OFI system. Fig. 3 shows the theoretical relationship between the amplitude of a plane acoustic wave in Pascals that passes through the laser beam, and the equivalent displacement in nanometers, as calculated using the model in [START_REF] Ciddor | Refractive index of air: new equations for the visible and near infrared[END_REF]. Plotted in a dashed blue line in the same figure we observe how the sensor's response the acoustic pressure follows the theoretical curve in a linear fashion.

Using [START_REF] Veng | Fringe Disappearance in Self-Mixing Interferometry Laser Sensors: Model and Application to the Absolute Distance Measurement Scheme[END_REF] we calculate that a 11.26 Pa zero to peak acoustic pressure wave corresponds to an equivalent displacement of 1.64 nm zero to peak. The lowest detectable amplitude in our experiment is 0.4 Pa zero to peak and is calculated to be 58 pm in equivalent displacement zero to peak. For each acoustic power we use the accelerometers to estimate relative displacement between the mirror and the laser. We observe a modulation of the relative distance at the acoustic frequency which is 3 kHz. In Fig. 3 we observe how the relative amplitude between the laser mount and the mirror mount descends linearly with the acoustic pressure. When injecting the 11.26 Pa zero to peak pressure wave into the system, we observe a 0.74 nm zero to peak relative movement between the laser mount and the mirror mount. As such we consider the vibrational levels too strong to be negligible. The average relative movement of the system is 0.038 nm zero to peak at the lowest detectable acoustic power of 0.4 Pa zero to peak. Therefore we calculate a new lower limit of detection for the system by adding the measured vibrational displacement to the theoretical equivalent displacement, and we obtain 96 pm zero to peak, that is similar to the lowest reported values for vibration OFI sensing in the literature. Thus we account for the possible perturbation due to vibrations in the system.

IV. CONCLUSION

We have studied the response of an OFI sensing system measuring the acousto-optic effect for different amplitudes of acoustic pressure variation. A model was established to convert the measured acoustic pressure amplitude to an equivalent displacement. Using planar acoustic waves made it possible to control and accurately measure the pressure p along the optical axis of the laser. The system was designed so that the acoustic waves were completely planar in order to calculate the equivalent displacement using the model.

We have found that there is a linear relationship between the optical power variation and the acoustic pressure variations. Due to vibrations in the system the relative distance between the laser mount and the reflector mount is modulated at the acoustic frequency. We have determined the lowest detectable acoustic amplitude detectable, taking into account these perturbations. The lowest detectable optical path variation we are able to detect with our system is 96 pm when the acoustic wave pressure variation is 0.4 Pa. The relative movement due to vibrations is 38 pm and the calculated equivalent displacement is 58 pm. In [START_REF] Giuliani | Self-mixing laser diode vibrometer[END_REF] the lower threshold of vibration reported is 0.1 nm, giving our system essentially the same performances. In order to eliminate the effect of the vibrations the experiment is to be improved so that the laser and mirror mounts no longer are in contact with the waveguide. The acoustic generator should be placed on vibration absorbent foam in order to eliminate vibrations propagating through the optical table.

A future perspective is to continue the development of our model to make it compatible with spherical free-field acoustic waves. Such a model would be of great advantage when studying the acousto-optic effect using OFI, and could even be applied to previous studies already published. In doing so we could, using OFI, create a non-contact, cost effective and accurate measurement instrument for the quantification of freefield acoustic pressure waves. As the instrument would be noncontact, the wave would not be perturbed as it is measured.

Fig. 1 :

 1 Fig. 1: Top view of experimental setup.

Fig. 2 :

 2 Fig. 2: Blue line: Sensor signal power for acoustic pressures. Dashed black line: Noise level of the sensor.

Fig. 3 :

 3 Fig. 3: Black solid line (left axis): Simulated equivalent displacement for given pressure variations. Black striped line (left axis): Differential movement between laser mount and mirror mount for given pressure variations. Blue striped line (Right axis): Sensor response for measured acoustic pressure

V. ACKNOWLEDGMENTS This work has been performed with the LabCom CapIROa joint laboratory with the Toulouse INP, The LAAS -CNRS and ACOEM. Grants and funding provided by the Agence Nationale de la Recherche (ANR) the Association Nationale Recherche Technologie (ANRT).