Jean-Pierre Da Costa

Frédéric Galland

Antoine Roueff

Christian Germain

Antoine Roue

Unsupervised segmentation based on Von Mises circular distributions for orientation estimation in textured images

Keywords: Texture, orientation, circular data, Von Mises, segmentation, polygonal active grid, minimum description length

come

Unsupervised segmentation based on Von Mises circular distributions for orientation estimation in textured images

INTRODUCTION

Texture is a keystone of many image analysis frameworks. Whether for classication, for indexation or for segmentation applications, the importance of texture descriptors cannot be overlooked when dealing with natural images. Among such textural descriptors, pattern orientation has proved to be crucial for human image perception.

1, 2 Besides, machine vision systems have also beneted from orientation and anisotropy-based features in many domains (e.g.

).

Contact: jean-pierre.dacosta@ims-bordeaux.fr, frederic.galland@fresnel.fr 1 In particular, some papers have reported the use of pattern orientations and texture anisotropies for segmentation purposes. To mention a few, Bigün and Du Buf 6 used complex moments in Gabor space to capture the symmetries of the local image content within a quad tree segmentation approach. Rousson et al. [START_REF] Rousson | Active unsupervised texture segmentation on a diusion based feature space[END_REF] and Luis Gar- cía et al., 8 inspired by tensor-valued image analysis techniques, used the structure tensor 9, 10 as a textural feature driving level-sets segmentation frameworks. Recently, the same authors proposed a combined segmentation approach using both tensor and color information within a common minimization framework. [START_REF] De Luis-García | Texture and color segmentation based on the combined use of the structure tensor and the image components[END_REF] In prior works, Da Costa et al. [START_REF] Da Costa | Segmentation of high resolution remote sensing images : application to the automatic delineation of vine elds[END_REF][START_REF] Da Costa | Delineation of vine parcels by segmentation of high resolution remote sensed images[END_REF] have also taken advantage of texture anisotropy for the segmentation of remote-sensed images.

Anisotropy was captured either by spectral Gabor features or by spatial anisotropy descriptors, and was used either in a region growing segmentation framework [START_REF] Da Costa | Segmentation of high resolution remote sensing images : application to the automatic delineation of vine elds[END_REF] or simply in a thresholding/regularizing scheme. [START_REF] Da Costa | Delineation of vine parcels by segmentation of high resolution remote sensed images[END_REF] In the context of manuscript analysis, Grana et al. [START_REF] Grana | Describing texture directions with von mises distributions[END_REF] proposed a classication algorithm to discriminate between sub blocks of text, illustrations and background using anisotropy. Orientation distributions inside sub blocks were captured by grey level autocorrelation and modelled using mixtures of Von Mises density functions.

15

In this paper, segmentation is not considered as an aim of image analysis but as a means of estimating orientations in textures. Our contribution is a segmentation-based parametric approach for the estimation of orientations in directional textures, i.e. textures exclusively composed of oriented patterns. [START_REF] Le Pouliquen | A new adaptive framework for unbiased orientation estimation[END_REF][START_REF] Michelet | Estimating local multiple orientations[END_REF] The approach relies on a rst step of local estimation of the orientations in the image using derivative techniques. [START_REF] Le Pouliquen | A new adaptive framework for unbiased orientation estimation[END_REF] The obtained orientation eld image is segmented into homogeneous regions, using the Minimum Description Length (MDL) segmentation technique formerly presented in 18 and generalized in this paper to π-periodic data modeled with circular Von Mises distributions. 15 A nal orientation estimation is performed by averaging low level orientations inside the segmented regions.

The paper is organized as follows. Section 2 is a methodological presentation dealing with orientation estimation and circular data modelling. Section 3 presents a generalization of the MDL segmentation technique to Von Mises circular data. Section 4 reports an experimental evaluation of the segmentation-based framework for orientation estimation on synthetic and real textured images. The last section provides some conclusions.

ORIENTATION AND TEXTURES

Textures of interest

As orientation relates to visual perception, dening the concept of image orientation is not trivial. The notion of orientation is often closely linked to the kind of patterns observed in the images. In this paper, only the case where a unique orientation (i.e. the orientation of the underlying oriented pattern) is present at any pixel of the image will be considered. This corresponds to the linearly symmetric model, 10 to the {i1D} case i.e. intrinsic dimensionality is 1 according to Krieger and Zetzsche's nomenclature 19 or to the mono-directional texture case mentioned by Michelet et al. [START_REF] Michelet | Estimating local multiple orientations[END_REF] The limitation to unique orientations excludes some textures where superimpositions or occlusions occur, 20 but such textures are beyond the scope of this paper.

Orientation Estimation

It is commonly accepted that anisotropy perception is closely linked to the scale at which an image is considered.

The multi-scale nature of orientation information has therefore been taken into account in many vision systems both in spectral approaches 6, 10 and in spatial approaches 5 which usually attempt to capture anisotropy properties at various scales. However, in this paper we chose to focus on local orientations only, thus restricting orientation estimation to local pixel neighborhoods. The objective is twofold. First, keeping orientation estimates local limits border eects at the boundaries of regions with dierent orientations. Secondly, as local estimates can be processed by spatial lters operating within very small pixel neighborhoods, strong dependences between orientation estimates on neighboring pixels are avoided.

In the case of mono-directional textures, the use of local derivatives 21, 22 is a natural approach for orientation estimation. In this paper we use the formerly proposed 16 gradient operator GOP. Similar to standard 3 × 3 operators (e.g. Sobel lter), GOP is part of a panel of lters designed to limit both bias and noise sensitivity while keeping the estimation process as local as possible. Local texture orientation is orthogonal to the orientation of the gradient. It is dened within an interval of width π, for instance ] -π 2 , π 2 ], as two values that are π apart refer to the same orientation. Such orientation data are called axial data 15 i.e. π-periodic circular data.

Nevertheless, in the case of noisy images, orientation estimations are very scattered and cannot be used directly. For example, the estimation of orientations inside the textured image of g. 1.a is shown on g. 1.b, using the GOP operator. In order to reduce the noise in these images, a classical solution consists in averaging the orientations in a small neighborhood V (x, y) around (x, y). In this paper, V (x, y) is dened as a sliding square window containing N V = d × d pixels and centered on pixel (x, y). Let θ(x, y) be the orientation estimated at pixel (x, y). Thus the double angle α(x, y) = 2θ(x, y) is a 2π-periodic variable. Its averaging α can be computed inside V (x, y) as suggested by Mardia and Jupp: 15

2θ(x, y) = α(x, y) = tan -1 S(x, y) C(x, y) (1)
where

   C(x, y) = 1 N V (u,v)∈V (x,y) cos α(u, v) S(x, y) = 1 N V (u,v)∈V (x,y) sin α(u, v) (2)
with N V the number of pixels in V (x, y). In equations 1-2, angular values are replaced by unit vectors (cos α(u, v), sin α(u, v)) which are averaged within V (x, y). This avoids erroneous estimations provided by classical averaging when angular data close to the limits of the orientation interval occur (-π or π for instance). The resulting orientation provided by this averaging technique (denoted AGOP in the following for Averaged GOP operator) is closely related to the orientation of the rst eigen vector of the structure tensor 9,10 which is another classical approach to get a regularized orientation estimation.

The averaging θ(x, y) = α(x, y)/2 of the orientation image of g. 1.b obtained using squared sliding windows V (x, y) with dierent side sizes d (and thus N V = d × d pixels) are reported on g. 1.c (d=7) and g. 1.d (d=11).

It is notably shown on the left part of g. 3 that the Mean Average Deviation (MAD) between the averaged estimated orientation {θ(x, y)} (x,y) and the true orientation {θ GT (x, y)} (x,y) (GT for Ground Truth) used to generate the data (see g. 1.h), is minimal for d=7, using the following denition for the MAD:

M AD = 1 N (x,y)∈Image min(|θ GT (x, y) -θ(x, y)|, π -|θ GT (x, y) -θ(x, y)|) ( 3
)
where N is the number of pixels in the image. Of course, the value d which minimizes the MAD, denoted d M AD in the following, highly depends on the image considered and on the noise level. For instance, d M AD =11 when the noise level is increased by a factor 2 (σ = 100) as shown in g. 2.

In real applications, where the ground truth is unknown and thus quality criteria such as the MAD cannot be calculated, the estimation of the size of this sliding window is a dicult problem. In this paper, instead of trying to nd the best window size, it is proposed to directly recover the dierent homogeneous regions in the image by decomposing the orientation eld into a patchwork of regions with homogeneous orientations. Inside each homogeneous region, the whole pixels can then be used to rene the orientation estimation. The estimation of the orientation is thus averaged on the biggest homogeneous sample of pixels, but without being corrupted by the presence of inhomogeneities near borders.

PARTITION INTO REGIONS WITH HOMOGENEOUS ORIENTATIONS

Let θ = {θ(x, y), (x, y) ∈ [1, N x ] × [1, N y ]} be an image of orientations composed of N = N x × N y pixels. This image is assumed to be composed of a patchwork of R homogeneous regions Ω r (r ∈ [1, R]).
In other words, the pixel orientations in each region Ω r are assumed to be independent realizations of a probability density function (PDF) P [βr] , where β r denotes the parameter vector of the PDF. A partition is dened with the rule

(x, y) ∈ Ω r =⇒ θ(x, y) ∼ P [βr]
. The goal of image partitioning is then to retrieve the underlying partition w(x, y), so that w(x, y) = r ⇐⇒ (x, y) ∈ Ω r . The main diculty of such a task is that the PDF P In this section, the approach proposed in 18 (initialy devoted to the segmentation of Synthetic Aperture Radar (SAR) images) is generalized to the segmentation of orientation eld images. It relies on the minimization of a MDL criterion, [START_REF] Rissanen | Modeling by shortest data description[END_REF][START_REF] Rissanen | Stochastic Complexity in Statistical Inquiry[END_REF] and it has been demonstrated 18, 25 that in the case of Gaussian, Gamma, Poisson or Bernoulli PDFs, it can lead to ecient techniques to obtain a decomposition of the image into statistically homogeneous regions. In the present case of orientation images, such approaches have to be extended to circular data.

[βr] , r ∈ [1, R] are unknown.
3.1 The Von Mises PDF for circular data Among all the PDF proposed in the literature to model 2π-periodic circular data, 26 a lot of interest has been devoted to the Von Mises PDF, 15 which has also been used for image processing to model orientation data in textures. [START_REF] Grana | Describing texture directions with von mises distributions[END_REF][START_REF] Péron | Joint linear-circular stochastic models for texture classication[END_REF] This PDF is unimodal and symmetric and is entirely dened by two parameters: its mean orientation and its concentration. In this paper, the π-periodic PDF of the orientation θ(x, y) will be directly derived from this model, by considering that α(x, y) = 2θ(x, y) is distributed according to the standard 2πperiodic Von Mises PDF, leading to

P [µr,κr] (α) = 1 2πI 0 (κ r )
e κr cos(α-µr) .

(

) 4
where I 0 is the modied Bessel function of the rst kind and order 0, and where µ r ∈] -π, π] and κ r ≥ 0 are respectively the mean orientation and the concentration of α inside Ω r . The mean orientation of θ = α/2 is thus equal to µ r /2 (modulo π).

In the following, the purpose will thus be to generalize the MDL-based partitioning technique proposed in [START_REF] Galland | Minimum description length synthetic aperture radar image segmentation[END_REF][START_REF] Galland | Multicomponent image segmentation in homogeneous regions by stochastic complexity minimization[END_REF] to such circular data.

Determination of the MDL criterion adapted to Von Mises data

According to 18, 25 the MDL criterion ∆(θ, w), which represents the code length needed to encode the whole image using the partition w, is the sum of three terms:

∆(θ, w) = ∆ G (w) + ∆ P (β|w) + ∆ L (θ|w, β) (5)
where ∆ G is the code length needed to encode the partition w, ∆ P is the code length needed to encode the PDF parameters β r in each region Ω r dened by w, and ∆ L is the code length needed to encode the pixel orientation

θ(x, y) inside each region Ω r knowing β = {β r } r∈[1,R] .
Let us now detail the expression of these three terms, as a function of the partition w.

Following 18 it is proposed to model the partition w with a polygonal grid, i.e. a set of nodes linked by segments to dene the dierent regions Ω r . In this case, the code length ∆ G needed to encode a given partition w can be approximated with 18 ∆ G (w) = n(log N + log p) + log p + 2p + p log(4 dx dy )

where p is the number of segments of the polygonal grid, n the minimum number of Eulerian graphs needed to describe the grid, and dx and dy the mean lengths of the horizontal and vertical projections of the grid segments and N the number of pixels of the image (see 18 and 28 for details).

The code length needed to encode the PDF parameters β r inside each region can be approximated with: 24 ∆ P (β|w

) = R r=1 card(β r ) 2 log(N r ) (7)
where N r is the number of pixels inside region Ω r and where card(β r ) is the size of the parameter vector β r , i.e.

the number of PDF parameters that need to be encoded. In the present case, the parameter vector of the Von Mises PDF is β r = (µ r , κ r ), and thus card(β r ) = 2.

The last term ∆ L represents the code length needed to describe the pixel orientation θ(x, y) inside each region using an entropic code:

29

∆ L (θ|w, β) = R r=1 -L(Ω r |β r ) (8)
where L(Ω r |β r ) is the log-likelihood of the sample {α(x, y), (x, y) ∈ Ω r }, knowing the PDF parameter β r . In the case of Von Mises PDF, one directly obtains: 15

L(Ω r |β r ) = N r log 2π + κ r R r cos(α r -µ r ) -log I 0 (κ r ) (9)
where

α r = tan -1 Sr Cr R r = C r 2 + S r 2 (10)
and

   C r = 1 N (x,y)∈Ωr cos α(x, y) S r = 1 N (x,y)∈Ωr sin α(x, y) (11)
Since β r = (µ r , κ r ) is generally unknown, it has to be estimated. Following 15 the parameters µ r and κ r of the Von Mises PDF can be estimated in the Maximum Likelihood (ML) sense, leading to the following ML estimates µ r and κ r :

µ r = α r κ r = A -1 (R r ) (12)
where A(κ) = I 1 (κ)/I 0 (κ) and I 1 (κ) = I 0 (κ). In the following, for numerical computation, the approximation

A -1 (R r ) ≈ (1.28 -0.53R r
2 ) tan(πR r /2) proposed in 30 (with maximum relative error of 0.032) will be used.

When replacing β r by its ML estimates β r = ( µ r , κ r ) in the log-likelihood of eq. 9, the following generalized log-likelihood L θ (Ω r | β r ) is obtained:

L(Ω r | β r ) = N r log 2π + κ r R r -log I 0 ( κ r ) (13)
The MDL criterion associated to a given partition thus contains neither unknown PDF parameter, nor tuning parameter to balance between its dierent terms. It is then required to nd the partition which minimizes this criterion. Following the optimization process proposed in, 18 starting from an initial regular thin grid delimiting 8×8 pixel regions (as shown in g. 6.c and 7.c), the algorithm consists in alternatively combining three grid modication steps in order to minimize the criterion: a region merging (RegM erg) step in which neighboring regions are merged (allowing the estimation of the number of regions), a node moving (N odM ov) step in order to estimate the contour location, and a node removal (N odRem) step in which nodes are progressively removed one by one (allowing to estimate the number of nodes and thus to regularize the contour). Typically, three cycles of RegM erg, N odM ov and N odRem steps are sucient to obtain satisfactory results (see 18 for details on the optimization process). For each modication of the grid tested by the optimization algorithm, the parameters µ r and κ r must be re-estimated in order to determine wether or not this modication leads to a decrease of the MDL criterion. Most of the time consuming of the algorithm is thus due to the computation of the 2D summations over the regions of C r and S r envolved in the calculation of µ r and κ r (see eq. 10, 11 and 12). In order to obtain a fast partitioning algorithm, the implementation proposed in 31 has been used, allowing one to replace these 2D

summations by 1D summations along the contours of the regions. Using such a technique, the partitioning of a 256×256 pixel image can then be obtained in less than one second on a standard personal computer (using C programming on a 3.2 GHz processor PC under Linux).

ORIENTATION ESTIMATION

First results on synthetic images

When applying this segmentation technique on the noisy orientation image of g. 1.b (obtained by applying a 3 × 3 GOP estimator), the segmentation results shown on gures 1.e,f are obtained. It is then proposed to directly estimate the orientation inside a region Ω r using the whole pixels of this region, i.e. using as large homogeneous samples as possible. The estimation of the orientation is then straightforward as the estimation of the mean orientation inside a region Ω r is θ r = α r /2. The orientation estimated inside each region obtained by the segmentation are shown on g. 1.g. It appears that these results are better than those obtained using sliding windows, not only in term of the precision of the contours, but also in term of quality since the MAD obtained The MAD values obtained using the segmentation-based estimation method have also been reported on these graphs.

using the proposed approach is lower than those obtained with sliding windows as shown on the left part of g.

3.

Comparison experiments

The method introduced in this paper was shown in the case of g. 1.a to outperform the AGOP gradient-based orientation estimation by local sliding window averaging. However, it still has not been compared to other standard methods for orientation estimation. In the right part of g. 3, we plot the evolution of MAD value as a function of scale for various orientation estimation methods. One of them is the steerable lter E2 32 with variable size. The second method is the IRON operator 17 whose size can also be tuned. These methods and the AGOP method discussed above are compared to the proposed segmentation-based method. The textured image used for this experiment is the synthetic image of g. 4a, which is made of several patches with dierent pattern frequencies and a noise level σ = 50.

In this case, orientation estimation based on segmentation leads to lower MAD than the other tested methods, whatever the scale, i.e. lter size, of the chosen operator. The corresponding results obtained with the segmentation based method and with each of the other tested methods (when selecting a posteriori the lter size minimizing the MAD) are presented in g. 4.

However, at very high noise level, some frontiers between adjacent regions could be missed with the proposed segmentation technique, which can lead to a degradation of the precision of the estimated orientation in these regions. To illustrate this problem, the experiment of g. 4 has been reproduced on g. 5 with progressively increasing the standard deviation σ. Fig. 5 shows the evolution of the MAD as a function of σ for two estimation methods: the proposed segmentation based estimation and the best AGOP lter (i.e. the AGOP lter which either when no constraint on the region size is imposed (gray solid line) or when limiting the size of the regions that can be merged in the segmentation algorithm (black solid line).

minimizes the MAD). As soon as σ > 90, the performances of the proposed segmentation based estimation technique rapidly decreases and are outperformed by the best AGOP lter.

Nevertheless, the robustness of this segmentation based technique can be increased at high noise level by limiting the size of the regions in the segmentation step. Indeed, as soon as the regions contain enough pixels to obtain an accurate estimation of the orientation inside the region, increasing the size of the region is not useful anymore, notably since it can lead to the merging of poorly contrasted regions. Since the size of the AGOP lter has been limited to 21 × 21 pixels, it has been chosen to impose the same constraint in the segmentation algorithm. For that purpose, in the region merging step of the segmentation algorithm, the merging of two regions with sizes greater than 21 × 21 pixels (for both regions) is not allowed. The obtained results are reported on g. 5: adding such a constraint on the maximal size of the regions allows one to obtain equivalent performances as the segmentation based approach for noise level σ < 80, while keeping at least the performances of the AGOP technique for high noise levels.

Application

To illustrate the use of our approach, the algorithm is nally applied to the processing of microscopy images provided by Snecma Moteurs (Safran group) of a Titane base alloy shown in gures 6 and 7. The structures of interest in such images are the large central regions composed of needles of similar orientation. The objective is to compute needle orientation statistics inside these large regions while getting their contours. The results provided in g. 6 and 7 show that the large regions are relevantly segmented and that their orientation can be successfully retrieved.

CONCLUSION

In this paper, a new segmentation based parametric approach for the estimation of orientation in textures has been proposed. This approach relies on an unsupervised MDL based segmentation technique that has been generalized to the case of π-periodic circular data modeled with Von Mises pdf. Local orientations are computed on each image pixel using a derivative approach. Then pixels are mapped into regions by the segmentation algorithm. Pixels are then assigned the mean orientation of the region they belong to.

Through a series of experiments using synthetic data, the approach appears to process the data successfully in most cases. By maximizing the size of the regions where orientation averaging takes place, accurate orientation inside the regions are estimated and precise boundaries between regions are also determined. On the dierent images analyzed in this paper, this approach allows one to obtain lower MAD than with the standard orientation estimation techniques that have been tested. Among the various perspectives, a validation of this technique on The MAD values obtained using the segmentation-based estimation method have also been reported on these graphs.

Figure 1 .

 1 Figure 1. (a) Textured image obtained by adding to a non noisy image (with grey levels between 0 and 100) a zero-mean Gaussian noise with standard deviation σ=50. The orientation estimated using a 3 × 3 GOP estimator is shown on (b) and the smoothing of this orientation with a d × d sliding windows is shown for d = 7 (c) and d = 11 (d). The segmentation into homogeneous regions obtained with the proposed technique is shown on (e,f) (either displayed on (a) or (b)) and the associated estimated orientations are shown on (g). The orientation ground truth is shown on (h). In this case d M AD = 7.

Figure 2 .

 2 Figure 2. Same as gure 1, but with σ = 100. In this case d M AD = 11.

Figure 3 .

 3 Figure 3. Left: MAD value as a function of the square size d of the sliding windows, in the case of the texture of g. 1.a. Right: MAD value as a function of the scale parameter for various orientation operators applied to the texture of g. 4.a.

Figure 4 .Figure 5 .

 45 Figure 4. Comparison experiment: (a) Synthetic patchy image in the range 0 -100 with additive noise of standard deviation σ = 50. (b) Orientation estimated using a 3 × 3 GOP estimator. (c) Segmentation into homogeneous regions obtained with the proposed technique. (d) Orientation deduced from (c). (e) Best result obtained with the E2 steerable lter (size 11). (f) Best result obtained with the AGOP method (averaging on a 11 × 11 window). (g) Best result obtained with the IRON operator (size 15). (h) True orientation on each patch.

Figure 6 .

 6 Figure 6. (a) Microscopy image of a Titane base alloy (360 × 277 pixels) provided by Snecma Moteurs (Safran group). (b) Orientation estimated using a 3 × 3 GOP operator. (c) Initial grid delimiting 8 × 8 pixel regions used to initialize the segmentation algorithm. (d,e) Segmentation results obtained with the proposed approach (either displayed on (a) or (b)) -computational time: 1.4 s. (f) Corresponding orientations estimated inside each regions.

Figure 7 .

 7 Figure 7. Same as g. 6 but on another Titane base alloy microscopy image provided by Snecma Moteurs (Safran group).

Figure captions Figure 1 .

 captions1 Figure captions Figure 1. (a) Textured image obtained by adding to a non noisy image (with grey levels between 0 and 100) a zero-mean Gaussian noise with standard deviation σ=50. The orientation estimated using a 3 × 3 GOP estimator is shown on (b) and the smoothing of this orientation with a d × d sliding windows is shown for d = 7 (c) and d = 11 (d). The segmentation into homogeneous regions obtained with the proposed technique is shown on (e,f) (either displayed on (a) or (b)) and the associated estimated orientations are shown on (g). The orientation ground truth is shown on (h). In this case d M AD = 7.

Figure 2 .

 2 Figure 2. Same as gure 1, but with σ = 100. In this case d M AD = 11.

Figure 3 .

 3 Figure 3. Left: MAD value as a function of the square size d of the sliding windows, in the case of the texture of g. 1.a. Right: MAD value as a function of the scale parameter for various orientation operators applied to the texture of g. 4.a.

Figure 4 .

 4 Figure 4. Comparison experiment: (a) Synthetic patchy image in the range 0 -100 with additive noise of standard deviation σ = 50. (b) Orientation estimated using a 3 × 3 GOP estimator. (c) Segmentation into homogeneous regions obtained with the proposed technique. (d) Orientation deduced from (c). (e) Best result obtained with the E2 steerable lter (size 11). (f) Best result obtained with the AGOP method (averaging on a 11 × 11 window). (g) Best result obtained with the IRON operator (size 15). (h) True orientation on each patch.

Figure 5 .

 5 Figure 5. Variation of the MAD when the standard deviation σ of the noise added on the image of g. 4.a increases. The orientation is estimated using the best AGOP lter (dashed black line) and the proposed segmentation based approach, either when no constraint on the region size is imposed (gray solid line) or when limiting the size of the regions that can be merged in the segmentation algorithm (black solid line).

Figure 6 .

 6 Figure 6. (a) Microscopy image of a Titane base alloy (360 × 277 pixels) provided by Snecma Moteurs (Safran group). (b) Orientation estimated using a 3 × 3 GOP operator. (c) Initial grid delimiting 8 × 8 pixel regions used to initialize the segmentation algorithm. (d,e) Segmentation results obtained with the proposed approach (either displayed on (a) or (b)) -computational time: 1.4 s. (f) Corresponding orientations estimated inside each regions.

Figure 7 .

 7 Figure 7. Same as g. 6 but on another Titane base alloy microscopy image provided by Snecma Moteurs (Safran group).

ACKNOWLEDGMENTS

The authors would like to thank Snecma Moteurs (Safran group), for providing images and the GdR ISIS group of the French National Research Agency (CNRS) for nancial support.


N Roch

E Flurin

F Nguyen

P Morfin

P Campagne-Ibarcq

M H Devoret

B Huard

Widely tunable, non-degenerate three-wave mixing microwave device operating near the quantum limit

Keywords:

We present the first experimental realization of a widely frequency tunable, non-degenerate threewave mixing device for quantum signals at GHz frequency. It is based on a new superconducting building-block consisting of a ring of four Josephson junctions shunted by a cross of four linear inductances. The phase configuration of the ring remains unique over a wide range of magnetic fluxes threading the loop. It is thus possible to vary the inductance of the ring with flux while retaining a strong, dissipation-free, and noiseless non-linearity. The device has been operated in amplifier mode and its noise performance has been evaluated by using the noise spectrum emitted by a voltage biased tunnel junction at finite frequency as a test signal. The unprecedented accuracy with which the crossover between zero-point-fluctuations and shot noise has been measured provides an upper-bound for the noise and dissipation intrinsic to the device.

Three-wave mixing devices, i.e. non-linear circuits converting power among three microwave signals, are key elements of analog information processing in the microwave domain [START_REF] Pozar | Microwave engineering[END_REF]. However, they are based on dissipative components such as semiconductor diodes, or SIS tunnel junctions biased near the superconducting gap [START_REF] Tinkham | Introduction to superconductivity[END_REF]. The loss of signal limits their operation and also introduces noise above the minimum required by quantum mechanics [START_REF] Caves | [END_REF]4]. A non-degenerate mixing device with noise close to that minimum level was demonstrated recently [5,6]. However, the hysteresis preventing flux tunability for this 4-junction circuit severely limited possible applications to analog quantum signal processing. In this Letter, we show that by adding four inductances to the 4-junction loop, we can fully suppress the hysteresis and reach a 500MHz frequency tunability while operating close to the quantum limit. Our improvement of the device tunability by an order of magnitude is obtained without jeopardizing other advantages of non-degenerate 3-wave mixing.

An ideal non-degenerate three-wave mixing device in the microwave domain absorbs three signals at frequencies such that ω X + ω Y = ω Z with complex amplitudes A in X , A in Y and A in Z , respectively, and reemits signals at the same frequencies with amplitudes

A out X , A out Y and A out Z such that |A out X | 2 + |A out Y | 2 + |A out Z | 2 = A in X 2 + A in Y 2 + A in Z 2
, that is without internal dissipation. The device can operate in two power amplification modes: i) the photon gain mode, for which A in X , A in Z is the pump power providing the energy difference between photons at ω Z and photons at ω X . The Josephson Parametric Converter (JPC) [6], consisting of a ring of four Josephson junctions, can perform both functions. However, its operation has little tunability since the flux Φ ext applied through the ring has to be adjusted in the close vicinity of the value Φ 0 /2, where Φ 0 = h/2e is the flux quantum. In the present work, we consider a more general 3-wave mixing device in which 4 linear inductances are cross-linking the ring-modulator like the spokes of a wheel (see Fig. 1a). The hamiltonian of the ring is

A in Z 2 A in X 2 , A in
H = -1 2 E J sin (ϕ ext )ϕ X ϕ Y ϕ Z + 1 2 (E L /2 + E J cos ϕ ext ) ϕ X 2 + ϕ Y 2 + 1 2 (E L /4 + E J cos ϕ ext ) ϕ Z 2 + O(ϕ X,Y,Z 4 ) (1)
where the three spatial mode amplitudes ϕ X = ϕ 1 -ϕ 3 , ϕ Y = ϕ 2 -ϕ 4 and ϕ Z = ϕ 1 + ϕ 3 -ϕ 2 -ϕ 4 are gaugeinvariant, orthogonal linear combinations of the superconducting phases of the four nodes of the Josephson junction ring (Fig. 1b).

We will see below how these standing wave modes can be excited by the propagating mode amplitudes A in X , A in Y and A in Z and emit the amplitudes A out X , A out Y and A out Z . In the hamiltonian (1), E L = ϕ 0 2 /L is the energy associated with each of the inductances L, and

E J = ϕ 2 0 /L 0 J
is the Josephson energy of each tunnel junction. We also define the reduced flux quantum ϕ 0 = /2e and the dimensionless flux ϕ ext = Φ ext /4ϕ 0 threading each of the nominally identical 4 loops of the device. The first term of the hamiltonian is a pure 3-wave mixing term, while the two others are quadratic terms determining the effective inductance of modes X, Y and Z : E L /2 > E J , modes X and Y can be tuned by varying ϕ ext while retaining their stability: L -1 X,Y > 0 on the whole range of variation. However, there is a range of fluxes for which L -1 Z < 0 where the device departs from ϕ Z = 0 so that the expansion (1) is inappropriate. If the inductances are lowered even more such that E L /4 > E J , then all three modes of the device are stable for every value of ϕ ext , but at the expense of significant dilution of the non-linear term. In contrast, as E L is lowered below 2E J , dilution of non-linearity is minimized, but at the expense of the stability of the three modes. This is why the JPC, for which E L = 0, can operate only within a small range of values of ϕ ext forbidding any tunability of the device.

L -1 X,Y,Z = ϕ -2 0 ∂ 2 H/∂ϕ
We have tested this new, tunable, mixing element design, by inserting the ring into a resonant structure consisting of two λ/2 transmission line resonators coupled to the X and Y modes (Fig. 1d) as in Ref. [7] is non-resonant and excited through resonator X using a hybrid coupler (Fig. 1d). By varying the externally applied flux, it is possible to adjust the X and Y resonator frequencies given by

ω X,Y = ω 0 X,Y π 2 L λ/2 X,Y /2 π 2 L λ/2 X,Y /2 + L X,Y (ϕ ext ) , (2)
where ω 0 X,Y is the resonance frequency of the bare λ/2 resonator without a ring, L λ/2 X,Y = 2Z 0 /(πω 0 X,Y ) its lumped-element equivalent inductance [START_REF] Pozar | Microwave engineering[END_REF] and Z 0 its characteristic impedance. As long as E L /4 + E J cos ϕ ext > 0, the ring inductance L X,Y is given by

L X,Y (ϕ ext ) = ϕ 2 0 E L 2 + E J cos ϕ ext -1 . (3)
The device presented in Fig. 1c is realized in a single e-beam lithography step. The critical current of the Al/Al 2 O 3 /Al Josephson junctions was designed to be in the µA range. The wide geometric linear inductances cross-linking the ring are approximately given by µ 0 l where l = 100 µm is the length of each of the four meanders. According to theory, they should present negligible kinetic inductance [8]. The value of the ratio E L /E J = 3 ± 2 should favor the stability of the X and Y modes.

The device was operated in the photon gain mode. The phase and amplitude of the waves A out X and A out Y , relative to those of A in X and A in Y are measured with a vector network analyzer, for a whole set of pump tones A in Z . Turning off the pump tone first, we obtained the resonance frequency of both resonators as a function of flux (see Fig. 2a) as well as their half-maximum bandwidths B X = 39 MHz and B Y = 29 MHz. Unlike in the JPC, no hysteresis was found in the dependence of the resonance frequency on applied flux, confirming the stability of our device. However, two regimes must be distinguished in the data: that of the wide arches obeying (2) with a ring inductance given by (3) and that of the narrow arches for which E L /4 + E J cos ϕ ext < 0 and where the ring inductance depends precisely on the nonzero value of ϕ Z emerging from the broken symmetry along the Z mode. It is interesting to note that the two possible opposite values for ϕ Z in this regime give exactly the same resonance frequency. Besides, the fit of Fig. 2a does not take into account the perturbative effect of the parasitic inductances in series with the junctions. Using the full hamiltonian and these stray inductances, a complete agreement with the data can be obtained over the full flux variation range [9].

The power gain G of the device is defined as the ratio of the reflected power with pump on and off. The dependence of the gain on the pump power is shown on Fig. 3. Note in particular that a dynamical bandwidth B = 3.2 MHz is obtained for a gain of 20 dB. We checked that the parametric amplifier relation

√ G × B(G) = 2 B -1 X + B -1
Y -1 holds to less than a MHz of deviation for any pump power yielding a gain greater than 5 dB, for both signal and idler waves, as theory predicts [5].

As illustrated on Fig. 2b, the amplifier center frequency can be flux-tuned over 400 MHz which represents a range two orders of magnitude greater than the bandwidth at 20dB. Indeed, for each center frequency, we can find a reproducible set of applied flux, pump power and pump frequency yielding a gain higher than 20 dB and a dynamical bandwidth of B = 3 MHz (Fig 2b). No amplification was found in the domain of the narrow arches. While this observation cannot be explained directly by the expansion (1), it is consistent with the full hamiltonian that predicts the non-linear term to be significantly spoiled by spurious terms when ϕ Z = 0. The key point of our experiment is that we can still benefit, outside the range of the narrow arches, from a confortable tunable 3wave non-linearity. The tunability of this non-degenerate amplifier can therefore compete with the state of the art degenerate Josephson amplifiers [10][11][12][13][14][15] with the added benefits of pump-signal separation.

We now turn to dynamical range measurements which further characterize the non-linear operation of our device. For these measurements, we first calibrated the attenuation of the line named "Signal In" (Fig. 1d) with an accuracy of 3 dB[16]. We then measured the so-called 1 dB compression point of the amplifier mode of our device, which is the input power for which the gain is reduced by 1 dB. As presented on Fig. 2b, this maximal power ranges between -133 dBm and -118 dBm, corresponding to 3 and 95 photons per inverse dynamical bandwidth. The reduction in maximal allowed power occurs at lower frequencies where we have also observed that the pump power needed for a given gain is ∼ 30 dB lower than at higher frequencies. We believe that it could be explained by the pump frequency becoming, at lower signal frequencies, resonant with a mode of the crossed resonators. The device would hence depart from the stiff pump condition needed for parametric amplification with maximal dynamic range.

In a last series of experiments, the noise of our device was assessed by using the noise emitted by a voltage biased NIN tunnel junction as input signal. This noise, which is well-understood and therefore of predictable amplitude, plays the role of an in situ calibrated signal. At small electronic temperatures (k B T e ω S ), the noise from a tunnel junction presents two regimes as a function of voltage. For eV < ω S , zero-point fluctuations across the junction dominate with a power spectral density S p (ω S ) = ω S 2 , while for eV > ω S , electrons in the junction produce non-equilibrium shot noise and S p (ω S ) = eV 2 . The electronic temperature T e in the electrodes of the junction sets the sharpness of the crossover between these two regimes [17,18] as S p = S + p + S - p with

S ± p (ω) = 1 4 (eV ± ω) coth eV ± ω 2k B T e . (4)
Our experiment was performed using an aluminum junction kept in its normal state by permanent magnets close-by. We measured in situ a normal resistance of 5) and fitting an overall gain Gsys = 94.6 dB and an extra noise N add of 2.8 quanta coming both from the unavoidable quantum noise of the idler port (0.5 quanta) and the unwanted losses between the tunnel junction and the amplifier (2.3 quanta). The gain Gsys allows us to express this power spectral density in units of photon number or quantum.

43.9 Ω (measurement lines not shown on Fig. 1b). The output spectral density was recorded with a spectrum analyzer and averaged over a 2 MHz bandwidth around the center frequency of the amplifier (see Fig. 4a). Its dependence with bias voltage was obtained (Fig. 4b) for an amplifier gain of 23 dB with the same settings as in Fig. 3. The measured power spectral density is remarkably well described by an expression of the form:

S p m (ω S ) = G sys (S p + N add ω S ). (5)
In the shot noise regime, it is possible to calibrate the system gain G sys = dS p m /d(eV /2) = 94.6 dB from the NIN tunnel element to the spectrum analyzer including a possible attenuation from the element to the input port of the amplifier. Without any additional calibration, we extracted the apparent system added noise N add = 2.8 at the plateau (Fig. 4b). This number of quanta can be thought of as the standard half quantum attributable to the unavoidable quantum noise of the load at the idler port, and 2.3 quanta left which can be seen as an upper bound on the extra noise generated inside the device. On the other hand, an electronic temperature T e equal to the refrigerator mixing chamber temperature of 35 mK describes perfectly the crossover. It is worth emphasizing that the noise power of the total measurement setup is presented in Fig. 4 without any background subtraction and is therefore the full absolute system noise. In fact, there is a finite attenuation between the junction and the amplifier leading to an underestimation of the gain counted from the input of amplifier and hence to the ac-tual noise added by the device. Besides the unwanted insertion loss inherent to our type of low temperature measurement setup, the complex impedance of the junction itself is imperfectly matched [19]. Given the size of the junction ( 10 µm 2 ) and previous experiments on similar junctions, we estimated its capacitance to be in the 0.7pF -1 pF range. Using the resistance of the junction and the characteristic impedance of the amplifier, we calculated that the loss of signal due to the RC filtering of the junction noise leads to an apparent added noise between 1.3 and 2.1 photons. Our measurement thus improves the N add found by Bergeal et al. in that the measurement frequency and bandwidth are substantially higher [5]. It is straightforward to compare the noise measurement with and without our device. Turning off the pump tone, the same noise measurement using only a state-of-the-art HEMT amplifier at 4 K [START_REF] Weinreb | Microwave Symposium Digest[END_REF] yielded an apparent added noise 20 times larger than with the pump on. This translates into an acquisition time 400 times longer, keeping the same bandwidth.

In conclusion, we have shown that it is possible to realize with Josephson tunnel junctions a widely tunable, dissipation-less, non-degenerate 3-wave mixing element which processes microwave signals, adding a level of noise not significantly greater than the level of unavoidable quantum noise. Such an element could be useful in a certain number of analog quantum signal processing applications, like the feedback control of the state of a quantum bit [START_REF] Korotkov | [END_REF].

Y 2

 2 is the pump power providing the extra photon numbers in the re-emitted signals at frequencies ω X and ω Y , and ii) the pure up-conversion mode for which A in Y 2

Figure 1 :

 1 Figure 1: a. Device schematic: four linear inductances L cross-link a ring of four Josephson junctions. Each sub-loop is biased by a magnetic flux ϕextϕ0. For L < L 0 J /4, the current through the inductances is zero and the external flux phase biases the junctions to ϕext. b. The device is embedded at the intersection of four transmission lines and couples to spatial modes X, Y and Z represented as arrows. c. Optical microscope image of the ring modulator. The meanders in the center of the ring implement the four linear inductances from a. The stripes on the meanders are due to the fabrication process based on shadow evaporation. d. Simplified schematic of the setup used to characterize 3-wave mixing operation. The idler resonator (X) is excited through a 180 • hybrid coupler while the signal resonator (Y ) is single-ended. The noise emitted by the voltage biased tunnel junction in its normal state is amplified through the signal port.

  . The Z mode

Figure 2 :

 2 Figure 2: a. Dots: Measured resonance frequency ωY of the signal cavity as a function of flux applied to the ring modulator without pump. Solid line: fit of ωY using Eq. (2) with ω 0 Y /2π = 8.82 GHz, L = 49 pH, EJ = ϕ0 × 1.9 µA and including the known stray inductance around the loop 4LS = 200 pH (see ref.[9]). b. Reflection gain measured on the signal port as a function of frequency for various values of the flux indicated by the color lines in a. Pump parameters are optimized for each curve. The numbers on top represent the 1 dB compression point (maximum input power) expressed in input photon rate per dynamical bandwidth for six different working frequencies coded by color.

Figure 3 :

 3 Figure 3: Reflection gain of the phase-preserving Josephson amplifier observed on the signal Y (a) and the idler X (b) modes. The color bar indicates the pump power referred to the output of the generator. The pump frequency is ωZ /2π = 14.071 GHz and the flux is set to ϕext/2π ≈ 0.3125.

Figure 4 :

 4 Figure 4: a. Power spectral density as a function of frequency measured at the output of "Signal Out" for three settings: pump off and V = 0, pump on and V = 0 or V = 100 µV. The colored area represents the averaging range used in the right panel. b. Average power spectral density over a 2 MHz bandwidth around the center frequency of the amplifier as a function of bias voltage V . The solid line shows what is expected using Eq. (5) and fitting an overall gain Gsys = 94.6 dB and an extra noise N add of 2.8 quanta coming both from the unavoidable quantum noise of the idler port (0.5 quanta) and the unwanted losses between the tunnel junction and the amplifier (2.3 quanta). The gain Gsys allows us to express this power spectral density in units of photon number or quantum.

  2 X,Y,Z . The value ϕ ext =π/2 maximizes the strength of the mixing term. Provided that

Acknowledgments: Discussions with the Quantronics group at CEA Saclay, as well as with F. Mallet, F. Schackert and T. Kontos have been greatly useful. We gratefully acknowledge O. Andrieu and J.C. Dumont for technical support. The devices have been made within the consortium Salle Blanche Paris Centre. We thank D. Mailly for helping us dicing wafers at LPN-Marcoussis and the Quantronics group for metal evaporations. This work was supported by the EMERGENCES program Contract of Ville de Paris and by the ANR contract ULAMSIG.