Marcos William

Silva Oliveira

Núbia Rosa Da Silva

Antoine Manzanera

Odemir Martinez

Odemir Martinez Bruno

Feature Extraction on Local Jet Space for Texture Classification

Keywords: Local jet space, pattern recognition, texture, image analysis

come   L'archive ouverte pluridisciplinaire

Introduction

Local jets make a decomposition of the image, that is an intensity function, based on partial derivatives [START_REF] Florack | The Gaussian scale-space paradigm and the multiscale local jet[END_REF][START_REF] Manzanera | Local jet based similarity for NL-means filtering[END_REF][START_REF] Manzanera | Local jet feature space framework for image processing and representation[END_REF]. This approach can be seen as a model of receptive fields of human cortex [START_REF] Koenderink | Representation of local geometry in the visual system[END_REF]. Receptive fields are able to discriminate geometric characteristics of the visual scene like lines, edges, orientations, or curvatures. In [START_REF] Manzanera | Local jet feature space framework for image processing and representation[END_REF], a local jet space framework was developed to process and represent images extracting relevant information over the decomposition of the image by a set of functions, demonstrating its efficiency in optical flow estimation, Non-Local means filtering, background subtraction and image and object characterization. Moreover, this framework can be used in a wide range of applications, such as real time video processing, detecting, predicting and filtering.

In the context of image representation, the behavior of local jet can be used to highlight the geometric characteristics of texture patterns. Furthermore, the local jet decomposition creates a similarity space of image that can be explored to classify texture images. Texture features have been played a significant role in pattern recognition and image analysis, including applications using real and synthetic images in several fields, for example plant analysis and identification [START_REF] Pydipati | Identification of citrus disease using color texture features and discriminant analysis[END_REF][START_REF] Journaux | Plant leaf roughness analysis by texture classification with generalized Fourier descriptors in a dimensionality reduction context[END_REF][START_REF] Rossatto | Fractal analysis of leaf-texture properties as a tool for taxonomic and identification purposes: a case study with species from Neotropical Melastomataceae (Miconieae tribe)[END_REF][START_REF] Florindo | Brachiaria species identification using imaging techniques based on fractal descriptors[END_REF], medical image analysis [START_REF] Castellano | Texture analysis of medical images[END_REF][START_REF] Zhou | Polarization imaging for breast cancer diagnosis using texture analysis and SVM[END_REF][START_REF] Yu | Coregistered FDG PET/CTbased textural characterization of head and neck cancer for radiation treatment planning[END_REF] and many others. Aiming to increase the potential of traditional texture descriptors, our proposal is to enhance the texture feature extraction process using the local jet decomposition as initial transformation. The texture image is represented by a set of features composed by the result of the application of a texture descriptor on the local jet components.

The performance of the proposal is assessed over four databases: Brodatz [START_REF] Manjunath | Texture features for browsing and retrieval of image data, Pattern Analysis and Machine Intelligence[END_REF], Usptex [START_REF] Backes | Color texture analysis based on fractal descriptors[END_REF], Vistex [START_REF]Vision texture database, Maintained by the Vision and Modeling group at the MIT Media Lab[END_REF], and Outex [START_REF] Ojala | Multiresolution gray-scale and rotation invariant texture classification with local binary patterns[END_REF]; using four texture feature descriptors: Fourier, Gabor, Local Binary Pattern and Local Binary Pattern Variance. A principal component analysis (PCA) [START_REF] Jolliffe | Principal component analysis[END_REF] was performed to reduce the dimension of features and linear discriminant analysis (LDA) is used to perform the classification. In most experiments, the proposed approach obtained higher success rate compared to the feature descriptors applied without the local jet decomposition.

The remainder of this paper is organized as follows. Section 2 presents the theory of local jet. Section 3 describes how is the process of feature extraction in the local jet space as well as each of the texture extractors used in the experiments. Section 4 outlines the experiments to validate the proposal, detailing the databases and procedures for classification. The results and discussion of the classification experiments are reported in Section 5, and, finally, conclusions are presented in Section 6.

Local Jet

To extract the basic local representation of images, we have used the multiscale local jet, which can be used as a similarity space [START_REF] Manzanera | Local jet based similarity for NL-means filtering[END_REF], based on the approximation provided by the Taylor expansion. For a function f : R 2 → R, at order r, the Taylor expansion is [START_REF] Florack | The Gaussian scale-space paradigm and the multiscale local jet[END_REF] where c = (c 1 , c 2 ) is the residual. The notation of partial derivatives can be changed to

f (x + c 1 , y + c 2 ) = r k=0 k i=0 k i c k-i 1 c i 2 ∂ k f ∂x i ∂y k-i (x, y) + O( c r ),
f ij = ∂ k f ∂x i ∂y k-i . In image processing context, f : Ω ∈ R 2 →
R is the original image and the derivatives are estimated at a certain scale [START_REF] Florack | The Gaussian scale-space paradigm and the multiscale local jet[END_REF], by convolving the image with Gaussian derivatives, as

f σ ij = f ∂ i+j G σ ∂x i ∂y j , (2)
where G σ is the bi-dimensional Gaussian filter with standard deviation σ, which represents the estimation scale. Hence, the local jet space is the set {f σ ij |i+j r, σ ∈ S}. In this work, the scale parameter is set to 1 (i.e. only one scale is used), and the order of derivation is set to 2. Thus, the local jet space of image f : Ω ∈ R 2 → R is the collection {f 00 , f 10 , f 01 , f 20 , f 02 , f 11 }. Figure 1 shows this decomposition applied to a texture image.

Each image in local jet space highlights different geometric properties of scene image. Figure 2 shows the isolines of a texture image and the isolines of some local jet components. We can see that the peaks in original image is blurred in f 00 component in which flat regions can be better analyzed. On f 10 resulting image, the intensity of directional variation is evidenced and the intensities in second order components represent the measure of concavity or convexity of each pixel neighbor on the image surface, emphasizing the different local variations.

Texture Feature Extraction on Local Jet Space

The primary goal of the proposed approach was to emphasize the intrinsic content in texture images that can not be extracted directly from the original image. The local jet decomposition enhances geometric properties in texture images to improve the action of feature extractors.

The flow in Figure 3 summarizes the idea where local jet decomposition until the second order on Equation(2) is applied in all images of the database, i.e., six local jet components are calculated for each image. A texture feature extractor is used on these local jet components characterizing each image in six n-dimensional feature vectors. The dimension n of each feature vector is dependent on the feature extractor. The PCA transformation is applied to the resulting set of feature vectors of each texture descriptor. This procedure aims to find the directions of the feature space along which the variability of feature vectors is maximal, allowing to reduce the dimension of the feature space.

Different feature extractors are used in this approach based on their different analysis of geometric properties. The application of this methodology is not limited to these feature extractors. They were chosen due to the wellknown applications of these descriptors in the literature. In the remainder, each descriptor is presented.

Fourier descriptors

The Fourier descriptors was proposed originally by Cosgiff in 1960 [START_REF] Zahn | Fourier descriptors for plane closed curves[END_REF] to caracterize a periodic signal. For texture classification, Azencott et al. [START_REF] Azencott | Texture classification using windowed Fourier filters, Pattern Analysis and Machine Intelligence[END_REF] extract texture features based on windowed Fourier filters.

The basic idea behind Fourier descriptors is that texture features are extracted on the frequency domain. Then, this process is performed by the application of 2D Fourier transform on the texture image f (x, y) multiplied by (-1) x+y to shift the origin of F (u, v) spectrum to center of the area occupied by the frequency domain [START_REF] Gonzalez | Digital Image Processing[END_REF] as

F[f (x, y)(-1) x+y ] = F (u -M/2, v -N/2), (3)
where M ×N is the size of image. This procedure allows the use of concentric circles on (0, 0) to obtain the features.

Based on [START_REF] Gonzalez | Digital Image Processing[END_REF][START_REF] Journaux | Texture classification with generalized fourier descriptors in dimensionality reduction context: An overview exploration[END_REF], the feature vector is composed by two different characteristics extracted from spectrum: Ea is the energy extracted from sectors of circular rings, and Eb is the energy from circular regions, defined as

Ea ij = λ i+1 λ i θ j+1 θ j |F (λ, θ)| 2 dθdλ, (4)
Eb i = λ i 0 2π 0 |F (λ, θ)| 2 dθdλ, (5)
where (λ, θ) are the polar coordinates of points in the frequency domain, and {λ i } M i=1 and {θ j } N j=1 are the sets of parameters for the method with M radius and N orientations. The texture image is wholly represented by the Fourier descriptor vector

F D = [Ea 11 , ..., Ea M N , Eb 1 , ..., Eb M ]. (6)

Gabor filters descriptors

Gabor filters [START_REF] Manjunath | Texture features for browsing and retrieval of image data, Pattern Analysis and Machine Intelligence[END_REF][START_REF] Daugman | Gabor wavelets and statistical pattern recognition[END_REF][START_REF] Idrissa | Texture classification using Gabor filters[END_REF] may identify edges in images according to a set of orientations and scales, as well as local image microstructures. This last property is applied to characterize texture images.

Formally, a 2D Gabor filter is a bidimensional Gaussian function modulated by sine and cosine function [START_REF] Idrissa | Texture classification using Gabor filters[END_REF][START_REF] Daugman | Gabor wavelets and statistical pattern recognition[END_REF]. The Gabor functions g(x, y) and its Fourier transform G(u, v) can be defined as

g(x, y) = 1 2πσ x σ y exp - 1 2 x 2 σ 2 x + y 2 σ 2 y + 2πjW x , (7)
G(u, v) = exp - 1 2 (u -W ) 2 σ 2 w + v 2 σ 2 v , (8)
where σ u = πσx 2 and σ v = πσy 2 . Based on g(x, y), the filter set can be obtained through dilations and rotations, i.e.,

g mn (x, y) = a -m G(x , y ), with a > 1, m, n ∈ N, x = a -m (x cos θ + y sin θ), and y = a -m (-x sin θ + y cos θ), (9)
where θ = nπ N and N is the number of orientations. The feature extraction procedure is based on the convolution of the texture image by a bank of Gabor filters with M scales and N orientations previously defined [START_REF] Manjunath | Texture features for browsing and retrieval of image data, Pattern Analysis and Machine Intelligence[END_REF]. In the frequency domain, the spectrum energy E is a general descriptor. Then, the feature vector is composed by the spectrum energy of each convolution. The texture image is wholly represented by the Gabor filter descriptors vector

GF D = [E 11 , E 12 , ..., E 21 , ..., E M N ]. (10)

Local Binary Pattern

Proposed by Ojala et al. [START_REF] Ojala | A comparative study of texture measures with classification based on featured distributions[END_REF][START_REF] Ojala | Multiresolution gray-scale and rotation invariant texture classification with local binary patterns[END_REF], Local Binary Pattern (LBP ) is a simple yet efficient operator to analyze texture images. This operator characterizes the spatial structure based on local image binary patterns. For each pixel (i, j) in image, the LBP code is computed by

LBP P,R (i, j) = P -1 p=0 s(g p -g c )2 p , where s(x) = 1, if x ≥ 0, 0, if x < 0, ( 11
)
g c is the gray level of central pixel (i, j), g p is the gray level of the p neighbors (p = 0, ..., P -1) and R is the radius of the neighborhood. The matrix with the LBP code of each pixel is the LBP map of the image and the histogram of the LBP map globally represents the texture.

In [START_REF] Ojala | Multiresolution gray-scale and rotation invariant texture classification with local binary patterns[END_REF], Ojala et al. define a uniformity measure U

U (LBP P,R (i, j)) = |s(g p-1 -g c ) -s(g 0 -g c )| + P -1 p=1 |s(g p -g c ) -s(g p-1 -g c )|
(12) and a contrast and rotation invariant texture descriptor using

LBP riu2 P,R (i, j) = P -1 p=0 s(g p -g c ), if U (LBP P,R (i, j)) ≤ 2, P + 1, otherwise. (13)
In our study, we consider the histogram of original LBP code to describe the texture image, and the rotation invariance is left only to the descriptor LBP V as described below.

Local Binary Pattern Variance

Ojala et al. [START_REF] Ojala | Multiresolution gray-scale and rotation invariant texture classification with local binary patterns[END_REF] suggest the use of variance measure V AR P,R to incorporate the contrast of local image texture information to LBP riu2 P,R operator. Hence,

V AR P,R = 1 P P -1 p=0 (g p -µ) 2 , where µ = 1 P P -1 p=0 g p , (14)
is combined with LBP riu2 P,R to obtain a powerful rotation invariant operator for local image texture analysis.

However, Guo et al. [START_REF] Guo | Rotation invariant texture classification using LBP variance (LBPV) with global matching[END_REF] show limitations of using V AR P,R by having continuous values that need to be quantized. This situation leads to the issue of choosing the number of bins of the histogram. On the one hand, if too few bins are used, they could not be enough to represent the image. On the other hand, if too many bins are used, it could generate a sparse or unstable histogram. Alternatively, Guo et al. propose the use of V AR P,R as an adaptive weight in histogram calculation of LBP code. Let the texture image be a M × N matrix with maximal LBP pattern value equal K, the whole texture image is represented by the Local Binary Pattern Variance (LBP V ) histogram

LBP V P,R (k) = M i=1 N j=1 w(LBP riu2 P,R (i, j), k), k ∈ [0, K], (15)
w(LBP riu2 P,R (i, j), k) =

V AR P,R (i, j), if LBP riu2 P,R (i, j) = k, 0, otherwise. (16)

Experiments

In this section, we show the various experiments that were carried out to demonstrate the better performance of extracting features from local jet space for texture classification instead of extracting features directly from the original images. In our experiments, we have used four databases: Brodatz, Vistex, Usptex and Outex, described in Table 1 with number of classes, number of images in all databases, number of images per class and size of each image. The entirety of each database was used in the experiments. To demonstrate the efficiency of our proposed approach, we have extracted features from images without preprocessing and we have compared with the results of extracting features over the local jet space. Fourier descriptors, Gabor filters, Local Binary Pattern and Local Binary Pattern Variance, as explained before in Section 3, were used in the experiments. Moreover, another set of experiments were conducted to verify the robustness of the proposed approach in front of noise images. To achieve this goal, the Gaussian noise, which appears during the process of image acquisition caused due poor illumination, electronic circuit noise, among others. Three different versions of each database have been generated using intensities of the Gaussian noise with zero mean and variance σ equals to 0.01, 0.05, e 0.1. Figure 4 shows samples of the modified databases with noise. As the feature space is a subset of R 6n , a dimensionality reduction method called Principal Components Analysis (PCA) was used. PCA converts a set of features extracted from image in orthogonal components to each others, so there is no redundant information. The transformed vector is as large as the original vector, however, that vector is a set of linearly uncorrelated values, sorted according to the variance of the original data [START_REF] Jolliffe | Principal component analysis[END_REF]. In all cases, the classification was performed using Linear Discriminant Analysis [START_REF] Fisher | The use of multiple measurements in taxonomic problems[END_REF] over the feature space transformed by PCA. The validation scheme of classification was stratified 10-fold cross-validation [START_REF] Hastie | The Elements of Statistical Learning[END_REF]. The experiments of comparison (i.e. the traditional methodology of texture extraction) were also conducted with PCA transformation over the feature space and classification using LDA, under stratified 10-fold cross-validation scheme.

Results and Discussions

For pattern recognition purpose, one of the ways of evaluating a feature extractor is to look at the dimension of descriptors. The lower the dimension of features, the more efficient the feature extractor. For this reason we analyse the success rate according to the dimension of features for each feature extraction method used here compared to the proposed local jet approach combining feature extraction method. Furthermore, these plots are shown for each database used in the experiments. Thus, it will be shown 16 graphics, being four databases (Brodatz, Vistex, Usptex and Outex) and four feature extractors (Fourier descriptors, Gabor filters, Local Binary Pattern and Local Binary Pattern Variance), as can be seen at Figures 5678. It can also be observed in these figures that the behavior of the curves is similar with respect to the increase of the dimension of descriptors (x-axis). Moreover, the proposed approach has better success rates compared to traditional methodology with low dimension of descriptors.

The results can be summarized and analyzed in two modes: with a fixed dimension of descriptors or comparing the best results, regardless of the dimension of descriptors. Considering the same dimension, the proposed approach has success rate higher than the approaches used in traditional way, as detailed in Figures 5678and Table 2. Here, the best improvement is the application of local jet associated with the Fourier descriptor in Usptex database. If we consider only the best success rate, ignoring the dimension of descriptors, the success rate rises considerably according to Table 3. In this regard, the proposed approach was responsible for an average increase of 2.6 percent in success rate for all traditional methods.

It is noticeable that in almost all experiments using local jet space resulted in a higher success rate. Two main observations should be done, the first is related to results by combining the local jet space with LBP method, because this combination produces the best success rate for all cases. The second observation relates to the results obtained in the experiments with Usptex and Outex, that are very difficult databases, where the best result without local jet decomposition is 84.43% for Usptex and 82.82% for Outex and using our proposal it was obtained 91.52% and 86.73%, respectively, according to Table 2.

Another points to be highlighted concern the results of LBPV, Gabor and Fourier. According to results obtained when using LBPV, which does not show an effective improvement by using the local jet space, coming to have lower success rate in the case of Usptex database for low dimension of descriptors. This is due to LBPV exploits complementary information for local variance of spatial patterns. However, local jet breaches the geometric patterns in each component reducing the local variance of these patterns. In general, this behavior benefits feature descriptors which disjointly measure the texture patterns. On the other hand, when the texture descriptor analyzes the variation among the patterns this property is less relevant in the local jet space.

Despite the similarity in the geometric interpretation of Gabor filters and local jet, the results of the combination Local jet + Gabor show a considerable increase in the success rate for all tests, as shown at Table 3. This is due to the fact that when applying the local jet over the image, important directional characteristics are isolated on each component and Gabor descriptors extract this information more accurately. In the same direction, Fourier descriptors also make use of the isolated directional characteristics by collecting relevant information in the same frequency points in the Fourier spectrum.

In the experiments related to noise invariance, again, the methodology using local jet space resulted in a higher success rate in most of the cases, as shown at Table 4. The results presented in Table 4 consider the number of descriptors that produces the best success rate for each experiment. According to the increasing of the Gaussian noise variance, the traditional methods tend to decrease the success rate, mainly when using LBP and LBPV descriptors. On the other hand, at the same time when using the local jet space and the same feature descriptors, the decay of success rate is smoother than in the conventional approach, proving the robustness of local jet space dealing with noisy images.

Conclusions

In conclusion, we investigated the potential of local jet space to improve the texture characterization by the highlight of the geometric patterns obtained through partial derivatives. In order to carry out this task, we applied classical texture descriptors to the six local jet components, resulting in a feature space of R 6n . To validate the proposal, the experiments were performed in Brodatz, Vistex, Usptex and Outex texture databases using Fourier, Gabor, LBP and LBPV as texture descriptors. In addition, PCA was used for dimensionality reduction, and to compare the new approach with the texture descriptors without using local jets considering feature vectors with the same dimension.

The results from the comparison between texture features extracted from the local jet space and directly from the image indicate the importance of the former, resulting in higher success rates for texture classification. For the particular experiment with Fourier descriptors and Usptex database, we have found the greatest increase of success rate, from 74.67% to 87.28%, by using of local jet decomposition. More generally, the proposed approach was responsible for an average increase of 2.6% in success rate for all traditional methods when considering only the best success rate of each experiment. The combination of the local jet space and texture feature extractors, leads us to the results that are promising for the field of pattern recognition and image analysis.

Table 4: Success rates of texture classification for three databases modified from the original ones obtained through addition of Gaussian noise with zero mean and variance σ = 0.01, 0.05 e 0.1.

Figure 1 :

 1 Figure 1: Local jet decomposition of one Brodatz texture image based on Taylor expansion.

Figure 2 :

 2 Figure 2: Isolines of local jet decomposition. (a) A cropping of one image form Brodatz database and its original image isolines (b). Isolines of (a): (c) f 00 , (d) f 10 , (e) f 20 and (f) f 11 .

Figure 3 :

 3 Figure 3: Proposed method scheme. The local jet decomposition is applied to each image of the database resulting in six local jet components normalized in [0, 255]. Fourier, Gabor, Local Binary Pattern and Local Binary Pattern Variance were used to extract n features from each component generating a feature vector v i ∈ R 6n . The technique for dimension reduction PCA is applied to the set of vectors of dimension 6n for each texture descriptor to obtain the coordinates with greatest variability. The feature space is composed by these transformed feature vectors.

Figure 4 :

 4 Figure 4: Image samples from the modified databases with Gaussian noise. Each row shows one image from the Brodatz, Vistex, Usptex and Outex databases, from top to bottom. First column represents the original image and the remainder columns represent the images with Gaussian noise zero mean and variance σ = 0.01, 0.05 e 0.1, from left to the right.

Figure 5 :

 5 Figure 5: Dimension of descriptors × Success rate comparing the traditional features extraction and feature extraction on local jet space for Brodatz database.

Figure 6 :

 6 Figure 6: Dimension of descriptors × Success rate comparing the traditional features extraction and feature extraction on local jet space for Vistex database.

Figure 7 :

 7 Figure 7: Dimension of descriptors × Success rate comparing the traditional features extraction and feature extraction on local jet space for Usptex database.

Figure 8 :

 8 Figure 8: Dimension of descriptors × Success rate comparing the traditional features extraction and feature extraction on local jet space for Outex database.

Table 1 :

 1 Databases.

	Database	# Classes # Images # Images per class	# Image size
	Brodatz [25]	111	1110	10	200 × 200 pixels
	Vistex [14]	54	864	16	128 × 128 pixels
	Usptex [13]	191	2292	12	128 × 128 pixels
	Outex [26]	68	1360	20	128 × 128 pixels

Table 2 :

 2 Success rate of all experiments using the same dimension of descriptors. # represents the dimension of descriptors and * means 163 descriptors for Brodatz database.

	Success rate (%) / (± std)

Table 3 :

 3 The best success rate of each experiment. # represents the dimension of descriptors.

		Dimension of descriptors (#) / Success rate (%) / (± std)
	Method/Database	Brodatz	Vistex	Usptex	Outex
		# % (± std)	# % (± std)	# % (± std)	# % (± std)
	Fourier	23 93.69 (± 0.17)	24 94.56 (± 0.28)	18 85.01 (± 0.31)	25 85.32 (± 0.42)
	Local jet + Fourier 42 95.56 (± 0.23)	69 97.84 (± 0.14)	50 89.63 (± 0.24)	63 87.16 (± 0.43)
	Gabor	19 92.18 (± 0.35)	18 90.09 (± 0.25)	20 81.80 (± 0.36)	19 82.96 (± 0.62)
	Local jet + Gabor	30 93.31 (± 0.28)	30 95.28 (± 0.36)	30 87.81 (± 0.23)	29 85.46 (± 0.31)
	LBP	84 97.86 (± 0.16)	60 98.87 (± 0.17)	80 90.93 (± 0.20)	73 85.51 (± 0.41)
	Local jet + LBP	216 99.02 (± 0.19) 91 99.86 (± 0.16) 104 94.29 (± 0.23) 56 88.98 (± 0.22)
	LBPV	10 87.99 (± 0.34)	10 81.35 (± 0.46)	10 73.65 (± 0.34)	10 64.20 (± 0.43)
	Local jet + LBPV	15 90.32 (± 0.35)	15 89.12 (± 0.31)	15 74.65 (± 0.30)	15 63.06 (± 0.69)

Preprint submitted to Physica A: Statistical Mechanics and itsApplicationsMay 28, 2015


H B Silveira
email: hector.silveira@ufsc.br

P S Pereira Da

Silva P Rouchon
email: pierre.rouchon@mines-paristech.fr

P S Pereira Da Silva

P Rouchon

Stabilization of photon-number states via single-photon corrections: a first convergence analysis under an ideal set-up

This paper presents a first mathematical convergence analysis of a Fock states feedback stabilization scheme via single-photon corrections. This measurement-based feedback has been developed and experimentally tested in 2012 by the cavity quantum electrodynamics group of Serge Haroche and Jean-Michel Raimond. Here, we consider the infinitedimensional Markov model corresponding to the ideal setup where detection errors and feedback delays have been disregarded. In this ideal context, we show that any goal Fock state can be stabilized by a Lyapunov-based feedback for any initial quantum state belonging to the dense subset of finite rank density operators with support in a finite photon-number sub-space. Closed-loop simulations illustrate the performance of the feedback law.

I. INTRODUCTION

In [START_REF] Zhou | Field locked to Fock state by quantum feedback with single photon corrections[END_REF], a photon-number states (Fock state) feedback stabilization scheme via single-photon corrections was described and experimentally tested. Such control problem is relevant for quantum information applications [START_REF] Nielsen | Quantum Computation and Quantum Information[END_REF], [START_REF] Haroche | Exploring the Quantum: Atoms, Cavities and Photons[END_REF]. The quantum state ρ corresponds to the density operator of a microwave field stored inside a super-conducting cavity and described as a quantum harmonic oscillator. At each sample step k ∈ N, a probe atom is launched inside the cavity. The measurement outcome y k detected by a sensor is the energy-state of this probe atom after its interaction with the microwave field. Each probe atom is considered as a two-level system: either it is detected in the lowest energy state |g , or the highest energy state |e . Consequently, the measurement outcomes corresponds to a discrete-valued output y k with only two distinct possibilities: g or e. Similarly, the control inputs u k are also discrete-valued with 3 distinct possibilities: -1, 0, +1. The open-loop value u k = 0 corresponds to a dispersive atom/field interaction: it achieves in fact a Quantum Non-Demolition measurement of Fock states [START_REF] Brune | Manipulation of photons in a cavity by dispersive atom-field coupling: Quantum-nondemolition measurements and generation of "Schrödinger cat" states[END_REF]. The two other values u k = ±1 correspond to resonant atom/field interactions where the probe atom and the field exchange energy quanta: these values achieve single-photon corrections.

Although the feedback law proposed and implemented in [START_REF] Zhou | Field locked to Fock state by quantum feedback with single photon corrections[END_REF] considered imperfect detections on y k and delays in the control, here we focus on an ideal-set up, that is, detection errors and control delays have been disregarded. Theorem 2 shows that, by adding an arbitrarily small term to the Lyapunov function used in [START_REF] Zhou | Field locked to Fock state by quantum feedback with single photon corrections[END_REF], one ensures almost sure global stabilization of any goal Fock state for the closedloop ideal set-up. This is achieved by relying on an infinitedimensional Markov model of the ideal set-up that takes into account the back-action of the measurement outcome y k on the quantum state ρ k+1 .

Loosely speaking, in [START_REF] Zhou | Field locked to Fock state by quantum feedback with single photon corrections[END_REF], the control value u k at each sampling step k was chosen so as to minimize the conditional expectation of the Lyapunov function V (ρ k ) = Tr (d(N )ρ k ), where N is the photon-number operator, d(n) = (n -n) 2 and ρ = |n n| is the goal Fock state. However, in closed-loop, the difference between such V and its conditional expectation is not strictly positive: such V does not become a strict Lyapunov function in closedloop and additional arguments have to be considered to prove convergence. These additional arguments are related to Lasalle invariance. They are well established in a smooth context where the control u is a smooth function of the state ρ. This cannot be the case here since u is a discrete-valued control. In order to overcome such technical difficulties, we propose, similarly to [START_REF] Amini | Feedback stabilization of discrete-time quantum systems subject to non-demolition measurements with imperfections and delays[END_REF], to add the arbitrarily small term -ǫ ∞ n=0 ( n|ρ k |n ) 2 to V (ρ k ), where ǫ > 0. This slightly modified control-Lyapunov function becomes then a strict-Lyapunov function in closed-loop that simplifies notably the convergence analysis. Moreover, the developed convergence analysis is done in the infinite-dimensional setting in the following sense: we show that, for any initial density operator ρ 0 with a finite photon-number support (ρ 0 |n = 0 for n large enough), the closed-loop trajectory k → ρ k remains also with a finite photon-number support with a uniform bound on the maximum photon-number. This almost finitedimensional behavior simplifies the convergence analysis despite the fact that such condition on ρ 0 is met on a dense subset of density operators (Hilbert-Schmidt topology on the Banach space of Hilbert-Schmidt self-adjoint operators).

The paper is organized as follows. Section II presents the ideal Markov model of the experimental set-up of the controlled microwave super-conducting cavity reported in [START_REF] Zhou | Field locked to Fock state by quantum feedback with single photon corrections[END_REF] and precisely formulates the Fock state stabilization problem here treated (see Definition 1). Section III establishes the proposed solution to the control problem in two distinct parts. Firstly, Section III-A considers the case where the initial condition ρ 0 is a diagonal density operator (see Theorem 1).

Only the main ideas of the convergence proof are outlined. The technical details are given in Section V. Afterwards, in Section III-B, the main result of the paper is presented: the general solution is obtained from Theorem 1 for ρ 0 belonging to a dense subset (see Theorem 2). The simulation results are exhibited in Section IV. The proof of some intermediate results and computations required in Sections III and V are presented in Appendices B-G. Finally, the concluding remarks are given in Section VI.

II. IDEAL MARKOV MODEL

Denote by H the separable complex Hilbert space L 2 (C) with orthonormal basis {|n , n ∈ N} of Fock states (photonnumber). Hence, H = { n∈N ψ n |n , (ψ 0 , ψ 1 , . . . ) ∈ l 2 (C)}. Let D be the set of all density operators on H, that is, the set of trace-class, self-adjoint, non-negative operators on H with unit trace. The sample step, corresponding to a sampling period around 100µs, is indexed by k ∈ N = {0, 1, 2, . . . }, u k ∈ {-1, 0, 1} is the control, ρ k ∈ D the quantum state and y k ∈ {g, e} the measurement outcome. The ideal Markov model of the controlled microwave superconducting cavity used in [START_REF] Zhou | Field locked to Fock state by quantum feedback with single photon corrections[END_REF] is given by:

ρ k+1 =      ρ g k+1 = M g (u k )ρ k M † g (u k ) Tr(M g (u k )ρ k M † g (u k )) when y k = g, ρ e k+1 = M e(uk )ρ k M † e (u k ) Tr(M e(uk )ρ k M † e (u k )) when y k = e, (1)
where the measurements outcomes y k = g and y k = e occur with probabilities 1 p g,k = Tr M g (u k )ρ k M † g (u k ) and p e,k = Tr M e (u k )ρ k M † e (u k ) = 1 -p g,k , respectively, u k = 0 corresponds to a dispersive interaction of the launched atom with the cavity field (Quantum Non-Demolition measurement of photons)

M g (0) = cos φ0N +φR 2 , M e (0)= sin φ0N +φR 2 , (2)
when u k = +1 the atom enters the cavity in the state |e with a resonant interaction with the cavity field

M g (+1) = sin θ0 2 √ N √ N a † , M e (+1)= cos θ0 2 √ N + 1 , (3)
when u k = -1 it enters in |g with a resonant interaction

M g (-1) = cos θ0 2 √ N , M e (-1) = a sin θ0 2 √ N √ N , (4
) and φ 0 , φ R , θ 0 ∈ R are adjustable control parameters. For each u ∈ {-1, 0, 1}, M g (u) and M e (u) are (linear) operators on H defined in the obvious way 2 according to the definitions in Appendix A. They are indeed well-defined operators on H, despite the fact that a and a † are unbounded 1 As usual in quantum physics, it is here assumed that the measurement outcome y k = y cannot occur when Tr M y (u k )ρ k M † y (u k ) = 0, for y = g, e.

2 For instance, M g (+1)|n = sin(

θ 0 2 √ N )/ √ N √ n + 1|n + 1 = sin( θ 0 2 √ n + 1)|n + 1 .
In order for the definition of M e(-1) to be consistent, it is assumed sin(0)/0 = 1.

operators. It is clear that M g (u), M e (u) are bounded operators on H with M † g (u)M g (u) + M † e (u)M e (u) = I (identity operator), M e (-1) = M † g (+1) = a sin( θ0 2 √ N )/ √ N , and M g (-1), M g (0), M e (0), M e (+1) are self-adjoint. It is easy to see that if the initial condition ρ 0 is a density operator then, for all realizations of the ideal Markov process (1)-( 4), ρ k is a density operator for k ∈ N.

Notice that ρ = |n n| is a steady state of the Markov process (1)-( 4) with u k = 0, where n ∈ N is arbitrary. The control problem here treated is given as follows:

Definition 1: For the ideal Markov process (1)-( 4), the control problem is to find a feedback law u k = f (ρ k ) such that, given an initial condition ρ 0 and n ∈ N, the closed-loop trajectory ρ k converges almost surely towards the goal Fock state ρ = |n n| as k → ∞.

The almost sure convergence above is with respect to the probabilities amplitudes P n (ρ) = Tr (|n n|ρ) = n|ρ|n of ρ, that is, lim k→∞ P n (ρ k ) = P n (ρ) for each n ∈ N. In other words, lim k→∞ P n (ρ k ) = 1 and lim k→∞ P n (ρ k ) = 0 when n = n. The solution proposed in this paper for the control problem above is developed in the next section.

III. STABILIZATION OF FOCK STATES

Given any operator A: H → H, let A mn = m|A|n for m, n ∈ N. Hence, A nn is the n-th diagonal element of A, while A mn with m = n correspond to its "off-diagonal" elements. One says that the operator A is diagonal when A mn = 0 for all m, n ∈ N with m = n. One shall begin by solving the control problem given in Definition 1 in the particular case where the initial condition ρ 0 is diagonal (see Theorem 1 in Section III-A). Afterwards, in Section III-B, the solution to the general non-commutative case is presented (see Theorem 2): its solution relies essentially on the diagonal case.

A. Diagonal case

For each n * ∈ N, define3

D n * = {ρ ∈ D | ρ is diagonal and ρ|n = 0, ∀n > n * } . Consider the set D * = n * ∈N D n * ⊂ D. Note that D n * ⊂ D n * +1 , and that each element ρ of D * is "finite dimensional" in the following sense: ρ ∈ D is in D n * if and only if ρ = n *
n=0 ρ nn |n n|, and ρ ∈ D n * may be considered as an operator from H to the finite-dimensional space H n * = span{|0 , . . . , |n * }, or as a density matrix on H n * . One defines the functions n min : D * → N, n max : D * → N and n length : D * → N respectively by:

• n min (ρ) is the smallest n ∈ N such that ρ|n = 0; • n max (ρ) is the greatest n ∈ N such that ρ|n = 0; • n length (ρ) = n max (ρ) -n min (ρ). It is clear that, given ρ ∈ D * , one has ρ ∈ D n * if and only if n max (ρ) ≤ n * .
The next result exhibits the properties of the state ρ k of ( 1)-( 4) with respect to these functions.

Proposition 1: For every realization of the ideal Markov process ( 1)-( 4) with initial condition ρ 0 ∈ D * , one has that ρ k ∈ D * for all k ∈ N with:

• If u k = 0 or u k = -1, then n max (ρ k+1 ) ≤ n max (ρ k )
and n length (ρ k+1 ) ≤ n length (ρ k );

• If u k = +1, then n max (ρ k+1 ) ≤ n max (ρ k ) + 1 and n length (ρ k+1 ) ≤ n length (ρ k ).
Proof: See Appendix B. Take a goal photon-number n ∈ N. As in [START_REF] Amini | Feedback stabilization of discrete-time quantum systems subject to non-demolition measurements with imperfections and delays[END_REF], consider the following Lyapunov function V ǫ : D * → R defined as

V ǫ (ρ) = Tr (d(N )ρ) -ǫ n∈N ρ 2 nn , for ρ ∈ D * , (5)
where ǫ > 0 is a real number and d(n) = (n-n) 2 as defined in [START_REF] Zhou | Field locked to Fock state by quantum feedback with single photon corrections[END_REF]. The feedback law u: D * → {-1, 0, 1} is given by

u = f (ρ) Argmin υ∈{-1,0,1} E [V ǫ (ρ k+1 ) | ρ k = ρ, u k = υ] . (6)
Note that for each ρ ∈ D * and n * ≥ n max (ρ), d(N )ρ in ( 5) is a well-defined self-adjoint, non-negative, trace-class operator on H, by considering d(N ) as an operator on H n * and ρ as an operator from H to

H n * . Indeed, d(N )ρ = n * n=0 ρ nn (n -n) 2 |n n|. Thus, (5) is well-defined. More- over, since H n * is invariant under ρ ∈ D * for n * ≥ n max (ρ), it is clear that Tr (d(N )ρ) = Tr H n * (d(N )ρ)
, where on the right-hand side one considers ρ as an operator on the finitedimensional space H n * and the trace is taken over H n * .

We have the following convergence result when ρ 0 ∈ D * : Theorem 1: Let n ∈ N and ǫ > 0. In ( 2)-( 4), assume that φ 0 /π and (θ 0 /π) 2 are irrational numbers, and take φ R = π/2 -nφ 0 . Consider the closed-loop Markov process (1)-( 4) with u k = f (ρ k ), where the feedback law f is as in [START_REF] Nielsen | Quantum Computation and Quantum Information[END_REF]. Then, given any initial condition ρ 0 ∈ D * , one has that ρ k converges almost surely towards ρ = |n n| as k → ∞.

Its proof is decomposed into two steps: First Step. Choose n ∈ N and ǫ > 0. Let n 0 = n length (ρ 0 ), r 0 = n min (ρ 0 ). Then, there exists an integer m 0 > n 0 + r 0 + n + 1 (depending on n 0 , r 0 , n and ǫ) such that, for all closed-loop realizations ρ k , one has ρ k ∈ D m0 for k ∈ N. Second Step. Choose irrational numbers φ 0 /π and (θ 0 /π) 2 in ( 2)-( 4), and take φ R = π/2nφ 0 . In D m0 , V ǫ is a strict super-martingale: for all density operators ρ in D m0 , one has

E [V ǫ (ρ k+1 ) | ρ k = ρ, u k = f (ρ)]-V ǫ (ρ) = -Q Vǫ (ρ, f (ρ)),
where Q Vǫ (ρ, f (ρ)) ≥ 0, and Q Vǫ (ρ, f (ρ)) = 0 if and only if ρ = ρ. The almost sure convergence follows then from usual results on strict super-martingales for Markov processes with compact state spaces.

The complete proof of the two steps above is presented in Section V. The general case where the initial condition ρ 0 is not necessarily diagonal is treated in the next subsection.

B. General case

Consider, for each n * ∈ N,

D n * = {ρ ∈ D | ρ|n = 0, ∀n > n * } ⊂ D n * +1 ,
and let D * = n * ∈N D n * ⊃ D * . It is clear that ρ ∈ D is in D n * if and only if ρ = n *
m,n=0 ρ mn |m n|. Consequently, D * is a dense subset of D when D is endowed with the subspace topology induced from the Hilbert-Schmidt norm. Indeed, let J 2 be the complex Banach space of all Hilbert-Schmidt operators on H with the Hilbert-Schmidt norm B 2 = ( m,n∈N |B mn | 2 ) 1/2 , for B ∈ J 2 [START_REF] Reed | Methods of Modern Mathematical Physics: Functional Analysis[END_REF], [START_REF] Conway | A Course in Operator Theory[END_REF]. Since D ⊂ J 2 and ρ ∈ D n * has the form ρ = n * m,n=0 ρ mn |m n|, the density property of D * in D is clear.

One has that ρ ∈ D n * may be considered as an operator from H to the finite-dimensional space H n * , or as a density matrix on H n * . Hence, d(N )ρ is a well-defined trace-class operator on H, by considering d(N ) as an operator on H n * and ρ ∈ D n * as an operator from H to H n * . Indeed, d(N )ρ = n * m,n=0 ρ mn (m -n) 2 |m n|, and it is trace-class because its range is finite-dimensional [START_REF] Reed | Methods of Modern Mathematical Physics: Functional Analysis[END_REF], [START_REF] Conway | A Course in Operator Theory[END_REF]. Consequently, the Lyapunov function V ǫ in (5), the feedback in [START_REF] Nielsen | Quantum Computation and Quantum Information[END_REF]

• V ǫ (ρ) = V ǫ (∆ρ), for ǫ > 0; • α -1 M y (u)ρM † y (u) belongs to D * with ∆ α -1 M y (u)ρM † y (u) = α -1 M y (u)(∆ρ)M † y (u); • M y (u)(∆ρ)M † y (u) nn = M y (u)ρM † y (u)
nn , for all n ∈ N. In particular, α = Tr M y (u)(∆ρ)M † y (u) . Proof: See Appendix G. Now, let ǫ > 0 and ρ = |n n|, where n ∈ N. Assume that ρ 0 ∈ D * . Let ρ k , k ∈ N, be the corresponding closed-loop trajectory for a fixed realization of (1)-( 4) with feedback u k = f (ρ k ), where f is as in [START_REF] Nielsen | Quantum Computation and Quantum Information[END_REF]. It is immediate from the proposition above that:

• ρ k ∈ D * , for k ∈ N; • ∆ρ k ∈ D * , k ∈ N,
is the corresponding closed-loop trajectory of ( 1)-( 4) for the initial condition ∆ρ 0 , the same realization (and with the same transition probabilities p e,k and p g,k ), as well as the same feedback

u k = f (ρ k ) = f (∆ρ k ); • Tr (|n n|ρ k ) = Tr (|n n|∆ρ k ), for any n ∈ N.
From these arguments, Theorem 1 and the fact that ∆ρ = ρ, one immediately obtains the following generic solution to the control problem, that is, when the initial condition ρ 0 belongs to the dense subset D * of D:

Theorem 2: Let n ∈ N and ǫ > 0. In ( 2)-( 4), assume that φ 0 /π and (θ 0 /π) 2 are irrational numbers, and take φ R = π/2 -nφ 0 . Consider the closed-loop Markov process (1)-( 4) with u k = f (ρ k ), where the feedback law f is as in [START_REF] Nielsen | Quantum Computation and Quantum Information[END_REF]. Then, given any initial condition ρ 0 ∈ D * , one has that ρ k converges almost surely towards ρ = |n n| as k → ∞.

IV. SIMULATION RESULTS

This section presents the closed-loop simulation results concerning the application of Theorem 2 above to the ideal Markov process (1)-( 4). The quantum experimental results exhibited in [START_REF] Zhou | Field locked to Fock state by quantum feedback with single photon corrections[END_REF] used the following control parameter values in (2)-( 4): φ 0 /π = 0.252 and θ 0 /π ≈ 2/ √ n + 1. However, according to the assumptions in Theorem 2, φ 0 /π and (θ 0 /π) 2 should be irrational numbers. Hence, here one chooses φ 0 /3.14 = 0.252 and θ 0 /3.14 = 2/ √ n + 1. One takes ρ 0 = 15 n=0 |n n|/16 ∈ D * as the initial condition, n = 10 for the goal Fock state ρ = |n n|, and ǫ = 10 3 as the gain for the feedback u k = f (ρ k ) in ( 5)- [START_REF] Nielsen | Quantum Computation and Quantum Information[END_REF]. Figure 1 exhibits the simulation results for one closed-loop realization with such choices and a final sample step of 120. It shows: the dynamics of the populations of ρ k (top), the controls u k (middle) and the simulated outcomes y k (bottom). The populations of ρ k correspond to the following observables:

A 1 = n-1 n=0 |n n| (n < n), A 2 = |n n| (n = n), A 3 = n>n |n n| (n > n).
Therefore, one sees from the dynamics of the populations that ρ k converges to ρ as k → ∞, which is in accordance with Theorem 2. Note that n|ρ k |n ≈ 1 and u k = 0 for all k > 45.

Recall that Theorem 2 assumes that ǫ > 0. In order to further analyze the performance of the Lyapunov-based feedback law here proposed, we now make a comparison with the one used experimentally in [START_REF] Zhou | Field locked to Fock state by quantum feedback with single photon corrections[END_REF], which corresponds to take ǫ = 0 in (5), i.e. to disregard the term -ǫ n∈N ρ 2 nn . Figure 2 presents the simulation results for one closed-loop realization of such case. The control parameters, ρ 0 and n = 10 are the same as above. Note that n|ρ k |n ≈ 1 and u k = 0 for all k > 78. In order to make a comparison in terms of the speed of convergence, define the settling time k s to be the smallest k ∈ N such that n|ρ k |n > 0.9 for all k ≥ k. One has k s = 45 for the case ǫ = 10 3 above, and k s = 78 for ǫ = 0. Therefore, in the two realizations here simulated, the choice of ǫ = 10 3 reduced the settling time k s by 42% with respect to ǫ = 0. This behavior is typical on an average basis, thereby justifying the term -ǫ n∈N ρ 2 nn in (5). Table I shows the average value k s and the standard deviation σ of k s for ǫ ∈ {0, 0.1, 1, 10, 10 2 , 10 3 , 10 4 , 10 5 }, where a total of 5000 realizations were simulated for each ǫ. Notice that when ǫ is relatively large or relatively small in comparison to ǫ = 10 3 , the average settling time k s deteriorated. Furthermore, although for ǫ = 10 5 one has that k s increased by nearly 22% in comparison to ǫ = 10 3 , the standard deviation σ decreased by nearly 62%. Computer simulations have suggested that a choice of ǫ > 0 which may perhaps significantly improve k s generally depends on the initial condition ρ 0 and on the goal Fock state ρ = |n n|, and it has to be determined heuristically. respectively. Note that V ǫ = V + ǫW . Define:

V. PROOF OF THEOREM 1 (DIAGONAL CASE)

Proof of the First Step

: Let ǫ > 0. Define V : D * → R and W : D * → R as V (ρ) = Tr (d(N )ρ) , W (ρ) = - n∈N ρ 2 nn , (7)
• Q W (ρ, u) = W (ρ) -E [W (ρ k+1 ) | ρ k = ρ, u k = u], • Q V (ρ, u) = V (ρ) -E [V (ρ k+1 ) | ρ k = ρ, u k = u], • Q Vǫ (ρ, u) = V ǫ (ρ) -E [V ǫ (ρ k+1 ) | ρ k = ρ, u k = u],
for ρ ∈ D * and u ∈ {-1, 0, 1}. The proof of Theorem 1 is a straightforward consequence of the next proposition: Proposition 3: Let ǫ > 0 and n 0 , r 0 , n ∈ N. There exists an integer m 0 > n 0 + r 0 + n + 1 (depending on ǫ, n 0 , r 0 , n) such that, for each ρ ∈ D * with n length (ρ) ≤ n 0 , if n max (ρ) = m 0 , then

Q Vǫ (ρ, -1) > max {Q Vǫ (ρ, 0), Q Vǫ (ρ, +1)} . In fact, given ρ 0 ∈ D * , let n 0 = n length (ρ 0 ) and r 0 = n min (ρ 0 ). Note that n max (ρ 0 ) = n 0 + r 0 < m 0 . By Proposition 1, ρ k ∈ D * with n length (ρ k ) ≤ n 0 , for all k ∈ N. Since u = f (ρ) maximizes Q Vǫ (ρ, f (ρ)), Proposition 3
implies that when n max (ρ k ) = m 0 for some k ∈ N, then the input u k will be always be equal to -1, and hence Proposition 1 ensures that

n max (ρ k+1 ) ≤ n max (ρ k ) = m 0 . Therefore, n max (ρ k ) ≤ m 0 , k ∈ N, showing the First Step.
The following two lemmas are instrumental for showing Proposition 3. Their proofs are given in Appendix D and Appendix E, respectively.

Lemma 1: Given an arbitrary nonzero θ 0 ∈ R, fix any a ∈ R such that 0 < a < 1/2. For all nonzero N 0 , N ∈ N, there exists an integer N > N big enough such that,

0 < 1/2 -a ≤ sin 2 θ0 2 √ n ≤ 1/2 + a, for n = N , N + 1, . . . , N + N 0 -1. Lemma 2: Let ρ ∈ D * . Then: • |Q W (ρ, u)| ≤ 1, for each u ∈ {-1, 0, 1}; • Q V (ρ, 0) = 0; • Q V (ρ, +1) = - n∈N ρ nn [2(n -n) + 1] sin 2 θ0 2 √ n + 1 ; • Q V (ρ, -1) = n∈N ρ nn [2(n -n) -1] sin 2 θ0 2 √ n .
The proof of Proposition 3 is shown in the sequel. Proof: Let ǫ > 0 and n 0 , r 0 , n ∈ N. One has to show that there exists m 0 > n 0 + r 0 + n + 1 such that, if ρ ∈ D * with n length (ρ) ≤ n 0 , then u = -1 always maximizes Q Vǫ (ρ, u) whenever n max (ρ) = m 0 . From Lemma 2 and the fact that Q Vǫ = Q V + ǫQ W , to complete the proof it suffices to show that:

• If ρ ∈ D * is such that n length (ρ) ≤ n 0 and n max (ρ) ≥ n 0 + n, then Q V (ρ, +1) ≤ 0; • There exists m 0 > n 0 + r 0 + n + 1 such that Q V (ρ, -1) > 2ǫ, whenever ρ ∈ D * is such that n length (ρ) ≤ n 0 and n max (ρ) = m 0 . Note that Q V (ρ, +1) = - nmax(ρ) n=nmin(ρ) ρ nn [2(n -n) + 1] sin 2 θ0 2 √ n + 1 ,
for any ρ ∈ D * . Thus, if n length (ρ) ≤ n 0 and n max (ρ) ≥ n+n 0 , then n min (ρ) ≥ n, and hence the first claim is shown. Now, fix 0 < a < 1/2 and let 4 N ≥ 1 2 2ǫ 1/2-a + 2n + 1 . Applying Lemma 1 for N 0 = n 0 + r 0 + 1 and such choice of N , one gets N > N in which 0 < 1/2 -a ≤ sin 2 θ0 2 √ n , 4 As N is an integer, it follows that N ≥ n + 1.

for n = N , N + 1, . . . , N + n 0 + r 0 . Take m 0 = N + n 0 + r 0 . Let ρ ∈ D * with n length (ρ) ≤ n 0 and n max (ρ) = m 0 . Note that m 0 > n 0 + r 0 + n + 1 and n min (ρ) ≥ N + r 0 . From Lemma 2 and the inequality above for 1/2 -a, one obtains

Q V (ρ, -1) = m0 n=nmin(ρ) ρ nn [2(n -n) -1] sin 2 θ0 2 √ n ≥ m0 n=nmin(ρ) ρ nn [2(n -n) -1](1/2 -a) ≥ m0 n=nmin(ρ) ρ nn [2(N -n) -1](1/2 -a) = [2(N -n) -1](1/2 -a) m0 n=nmin(ρ) ρ nn .
Using the fact that m0 n=nmin(ρ) ρ nn = 1 and N >

Proof of the Second Step:

Let ǫ > 0. Recall that, by definition,

Q Vǫ = Q V + ǫQ W .
Using the same notation of the First Step, the central idea of the proof is to show that, given ρ ∈ D m0 , one has that Q Vǫ (ρ, f (ρ)) ≥ 0, and that Q Vǫ (ρ, f (ρ)) = 0 if and only if ρ = ρ. The following lemma is instrumental for the proof of such property. Its proof is presented in Appendix F.

Lemma 3: Assume that φ 0 /π is an irrational number in (2), and take φ R = π/2 -nφ 0 , where n ∈ N. Let ρ ∈ D * . Then:

• Q W (ρ, 0) ≥ 0, and Q W (ρ, 0) = 0 if and only if ρ = |n n| for some n ∈ N;

• Q W (ρ, +1) = Q W (ρ, -1) = 0 whenever ρ = |n n| for some n ∈ N. One has that m 0 > n, and (θ 0 /π) 2 is an irrational number by assumption. Recall that sin 2 (x) = 0 if and only if x = ℓπ, where ℓ is an integer. First we show that Q Vǫ (ρ, f (ρ)) = 0. By Lemma 2:

Q V (ρ, +1) = -sin 2 ( θ0 2 √ n + 1) < 0; Q V (ρ, -1) = -sin 2 ( θ0 2 √
n) < 0 when n > 0, and Q V (ρ, -1) = 0 when n = 0; and Q

V (ρ, 0) = 0. As Q W (ρ, u) = 0, for u ∈ {-1, 0, 1}, and u = f (ρ) maximizes Q Vǫ (ρ, u), one has that Q Vǫ (ρ, f (ρ)) = 0. Now, let ρ ∈ D m0 ⊂ D * . Since u = f (ρ) maximizes Q Vǫ (ρ, u), it follows that Q Vǫ (ρ, f (ρ)) ≥ Q V (ρ, 0) + ǫQ W (ρ, 0) = ǫQ W (ρ, 0) ≥ 0.
Suppose Q Vǫ (ρ, f (ρ)) = 0. Hence, Q W (ρ, 0) = 0, and so ρ = |n n| for some n ∈ {0, 1, . . . , m 0 }. It suffices to show that

Q Vǫ (ρ, f (ρ))| |n n| > 0 for n ∈ {0, 1, . . . , m 0 } with n = n. Assume that n > n. It is clear that u = f (|n n|) = -1 and Q Vǫ (ρ, f (ρ))| |n n| = [2(n-n)-1] sin 2 ( θ0 2 √ n) > 0. Assume now that n < n. Then, u = f (|n n|) = +1 and Q Vǫ (ρ, f (ρ))| |n n| = -[2(n -n) + 1] sin 2 ( θ0 2 √
n + 1) > 0. This completes the proof of the referred property.

The remaining part of the proof of the Second Step is a straightforward consequence of the standard stochastic convergence result below: Theorem 1,p. 195] Let Ω be a probability space and let W be a measurable space. Consider that X k : Ω → W , k ∈ N, is a Markov chain with respect to the natural filtration. Let Q: W → R and V : W → R be measurable non-negative functions with V (

Theorem 3: [5,
X k ) integrable for all k ∈ N. If E [V (X k+1 ) | X k ] -V (X k ) = -Q(X k ), for k ∈ N, then lim k→∞ Q(X k ) = 0 almost surely.
Indeed, let J 1 be the complex Banach space of all traceclass operators on H with the trace norm • 1 , that is, B 1 = Tr (|B|), where |B| √ B † B, for B ∈ J 1 . Recall that B ≤ B 1 and5 | Tr (AB) | ≤ A B 1 , for every B ∈ J 1 and each bounded operator A: H → H, where • is the usual operator norm (sup norm of bounded operators) [START_REF] Reed | Methods of Modern Mathematical Physics: Functional Analysis[END_REF], [START_REF] Conway | A Course in Operator Theory[END_REF]. Consider the subspace topology on D m0 with respect to J 1 . One has that the closed-loop trajectory ρ k , k ∈ N, is a Markov chain with phase space D m0 (with respect to the natural filtration and the Borel algebra on D m0 ). It is clear that D m0 is compact, and that Q ǫ and V ǫ -α ǫ are non-negative and continuous on D m0 , for all ǫ > 0, where α ǫ min ρ∈Dm 0 V ǫ (ρ). The theorem above implies that ρ k converges almost surely towards ρ as k → ∞ (with respect to the trace norm). This completes the proof of Theorem 1.

VI. CONCLUDING REMARKS

This paper provided a convergence analysis of Fock states stabilization via single-photon corrections under an ideal setup, that is, assuming perfect measurement detection and no control delays. In terms of convergence speed, the simulation results here presented have justified the inclusion of the term -ǫ n∈N ρ 2 nn in the Lyapunov-based feedback law (5)- [START_REF] Nielsen | Quantum Computation and Quantum Information[END_REF]. It is straightforward to verify that the convergence analysis developed in this paper remains valid for: (i) any other function d(n) in (5) satisfying d(n) = 0, d(n) is increasing for n > n and d(n) is decreasing for n < n; and (ii) ǫ > 0 dependent on n, that is, to take the term -n∈N ǫ n ρ 2 nn . However, it is an open problem how to choose the function d(n) and the gains ǫ n > 0 so as to achieve the best convergence speed.

Finally, the feedback law used in [START_REF] Zhou | Field locked to Fock state by quantum feedback with single photon corrections[END_REF], which corresponds to ǫ = 0, was tailored for an experimental set-up with measurement imperfections and control delays. The convergence analysis of such realistic situation will be investigated in the future.

a † |n = √ n + 1 |n + 1 . Note that these operators cannot be extended to H. Let f : N → R be a function. Define the operator f (N ): H n * → H n * by f (N )|n = f (n)|n , for each n = 0, . . . , n * . It is clear that f (N ) can be extented to H whenever f is a bounded function. Given f : N → R and an integer m, one defines g: N → R as: g(n) = f (n + m), when n + m ≥ 0; and g(n) = 0, when n + m < 0. One abuses notation letting f (N + m) stand for g(N ). Given two functions f, g: N → R, it is clear that f (N )g(N ) = g(N )f (N ) = (f g)(N ) and (f + g)(N ) = f (N ) + g(N ). Furthermore: aa † = N + I, a † a = N , af (N ) = f (N + 1)a, a † f (N ) = f (N -1)a † .

B. Proof of Proposition 1

2 :

 2 and n max can be extended to D * . Define the map ∆: D * → D * ⊂ D * as ∆ρ = nmax(ρ) n=0 ρ nn |n n|. Note that ∆ extracts the diagonal of ρ ∈ D * . It is easy to see that n max (∆ρ) = n max (ρ) and (∆ρ) nn = ρ nn , ρ ∈ D * . Moreover, ∆ρ = ρ when ρ ∈ D * . Other properties of the map ∆ are given in the next result: Proposition Let ρ ∈ D * , u ∈ {-1, 0, 1}, y = g, e. Take α = Tr M y (u)ρM † y (u) . Then: • Tr (Aρ) = Tr (A∆ρ), for every diagonal bounded operator A: H → H;

Fig. 1 .

 1 Fig. 1. Simulation of one closed-loop realization with gain ǫ = 10 3 : convergence of ρ k towards ρ (top), controls u k (middle), and outcomes y k (bottom). Notice that n|ρ k |n ≈ 1 and u k = 0 for all k > 45.

Fig. 2 .

 2 Fig. 2.Simulation of one closed-loop realization with gain ǫ = 0: convergence of ρ k towards ρ (top), controls u k (middle), and outcomes y k (bottom). Notice that n|ρ k |n ≈ 1 and u k = 0 for all k > 78. TABLE I AVERAGE SETTLING TIME ks AND STANDARD DEVIATION σ AS A FUNCTION OF THE GAINS ǫ, CONSIDERING 5000 REALIZATIONSǫ = 0 ǫ = 0.1 ǫ = 1 ǫ = 10 k s = 79.94 ks = 79.95 ks = 81.24 ks = 71.33 σ = 164.97 σ = 166.61 σ = 174.29 σ = 150.95 ǫ = 10 2 ǫ = 10 3 ǫ = 10 4 ǫ = 10 5 k s = 60.41 ks = 44.18 ks = 47.05 ks = 53.77 σ = 119.39 σ = 44.12 σ = 37.37 σ = 16.84

  a + 2n + 1 , one shows the second claim, thereby completing the proof of Proposition 3.

Fix 2 √ 2 √ n + 1 2 √ 2 √ 2 √ρnn cos 2 θ 0 2 √ n + 1 ρnn cos 2 θ 0 2 √ρn+1,n+1 sin 2 θ 0 2 √ n + 1

 22122221221 any ρ ∈ D * and let n ∈ N. In particular, ρ|n = ρ nn |n . It then follows from (2)-(4) that:M g (0)ρM † g (0)|n = ρ nn cos 2 φ0n+φR 2 |n , M e (0)ρM † e (0)|n = ρ nn sin 2 φ0n+φR 2 |n , M g (+1)ρM † g (+1)|n = 0, for n = 0, ρ n-1,n-1 sin 2 θ0 n |n , n ≥ 1, M e (+1)ρM † e (+1)|n = ρ nn cos 2 θ0 |n , M g (-1)ρM † g (-1)|n = ρ nn cos 2 θ0 n |n , M e (-1)ρM † e (-1)|n = ρ n+1,n+1 sin 2 θ0 M g (+1)ρM † g (+1) = nmax(ρ)+1 n=n min (ρ)+1ρn-1,n-1 sin 2 θ 0 n |n n|, (10)M e(+1)ρM † e (+1) = nmax(ρ) n=n min (ρ) |n n|,(11)M g (-1)ρM † g (-1) = nmax (ρ)n=n min (ρ) n |n n|, (12) |n n|.(13)By assumption, ρ 0 ∈ D * . Then, (1), (8)-(13) above and induction on k show the assertions in Proposition 1.C. Computation ofQ V (ρ, u)Fix any ρ ∈ D * and n ∈ N. Recall that V (ρ) = Tr (d(N )ρ), where d: N → R be given by d(n) = (n -n) 2 . Note that (1) implies that, for each u ∈ {-1, 0, 1},E [V (ρ k+1 ) | ρ k = ρ, u k = u] = Tr d(N )M g (u)ρM † g (u) + Tr d(N )M e (u)ρM † e (u) . (14)Take u = 0. From (8)-(9) in Appendix B, one hasE [V (ρ k+1 ) | ρ k = ρ, u k = 0] = Tr d(N )M g (0)ρM † g (0) + Tr d(N )M e (0)ρM † e (0) = Tr d(N ) M g (0)ρM † g (0) + M e (0)ρM † e (0) = Tr (d(N )ρ) = V (ρ).In particular, Q V (ρ, 0) = 0.(15)Now, take u = +1. Then, (14) above and (10)-(11) in

Note that if ρ = |n n| for some n ∈ N, then ρ ∈ Dn.

One also recalls that if A is a bounded operator on H and B ∈ J 1 , then AB, BA ∈ J 1 with Tr (AB) = Tr (BA).

VII. ACKNOWLEDGMENTS

The authors are indebted I. Dotsenko, M. Brune and J. M. Raimond for valuable discussions on the experimental feedback scheme.

was fully supported by CNPq (National Counsel for the Scientific and Technological Development), Ministry of Science and Technology, Brazil, P. S. Pereira da Silva was partially supported by CNPq. P. Rouchon was partially supported by Projet Blanc ANR-2011-BS01-017-01 EMAQS. H. B. Silveira is with Departamento de Automac ¸ão e Sistemas (DAS),

By summing and subtracting Tr sin

In particular,

Finally, take u = -1. Using (14) above and (12)-

By summing and subtracting Tr sin

In particular,

D. Proof of Lemma 1

Assume that N 0 is even (otherwise one may take N 0 + 1 instead of N 0 in this proof). Define the function η:

By definition, one has θ0

Using the definition of h and the symmetries 6 of the function sin 2 (•), it is easy to show that

Let ℓ ∈ N be even and big enough such that the following two conditions are simultaneously met:

Now, take N = ⌊η(ℓ)⌉-N0 2 +1 > N , where ⌊η⌉ denotes the greatest integer which is less or equal to η. By construction, η(ℓ) is in-between the points N + N 0 /2 -1 and N + N 0 /2, and hence it is in the interval

x. From the fact that φ ′ (x) = θ0 4 √ x , by the mean value theorem applied to the function φ and the second inequality in (20), one obtains

Then, the proof follows easily from (18), ( 19) and the fact that sin 2 (x -ℓπ/2) = sin 2 (x), for every even ℓ ∈ N.

E. Proof of Lemma 2

Proof of the first claim:

, where p g,k , p e,k ≥ 0 with p g,k + p e,k = 1.

is the difference of two numbers that are inbetween -1 and 0, one concludes that |Q W (ρ, u)| ≤ 1.

The second, third and fourth claims, are immediate from (15), ( 16) and (17) in Appendix C, respectively.

F. Proof of Lemma 3

Proof of the first claim:

, for y = g, e.

Hence, αρ g + (1 -α)ρ e = ρ, where α p g,k = Tr M g (0)ρM † g (0) . In particular, αρ g nn + (1 -α)ρ e nn = ρ nn , for n ∈ N. Note that, if α = 0, then M g (0)ρM † g (0) = 0, and so ρ e = ρ. Similarly, α = 1 implies ρ g = ρ. Thus, the identity αρ g nn + (1 -α)ρ e nn = ρ nn , for n ∈ N, still holds when α = 0 or α = 1. From (1), [START_REF] Reed | Methods of Modern Mathematical Physics: Functional Analysis[END_REF] and α = p g,k , one has

6 More precisely, sin

thereby showing the first part of the first claim. If ρ = |m m| for some m ∈ N with 0 < α < 1, then (8)-(9) in Appendix B imply that ρ g = ρ e = ρ, and so Q W (ρ, 0) = 0. Now, one shows that ρ = |m m| for some m ∈ N whenever Q W (ρ, 0) = 0. Suppose Q W (ρ, 0) = 0. Then, (21) implies that α = 0, or α = 1, or ρ g nn = ρ e nn for all n ∈ N with 0 < α < 1. Assume that α = 0. Hence, M g (0)ρM † g (0) = n∈N ρ nn cos 2 ( φ0n+φR

2

)|n n| = 0 by (8) in Appendix B. Suppose that ρ = |m m| for every m ∈ N. Thus, there exists n 1 , n 2 ∈ N with n 1 = n 2 , ρ n1,n1 > 0, ρ n2,n2 > 0. Recall that sin(x 1 ) = ± sin(x 2 ) if and only if

(assuming no division by 0). Now, using (1) and the above computations, one gets


Antoine Tran
email: antoine.[tran@ensta-paristech.fr

Antoine Manzanera
email: manzanera]@ensta-paristech.fr

A versatile object tracking algorithm combining Particle Filter and Generalised Hough Transform

Keywords: Object Tracking, Generalised Hough Transform, Particle Filter, Colour and Edge Features

This paper introduces a new object tracking method which combines two algorithms working in parallel, and based on low-level observations (colour and gradient orientation): the Generalised Hough Transform, using a pixel-based description, and the Particle Filter, using a global description. The object model is updated by combining information from a back-projection map computed from the Generalised Hough Transform, providing for every pixel the degree to which it may belong to the object, and from the Particle Filter, providing a probability density on the global object position. The proposed tracker makes the most of the two algorithms, in terms of robustness to appearance variation like scaling, rotation, non-rigid deformation or illumination changes.

I. INTRODUCTION

Object tracking is a challenging problem in computer vision, with applications in many areas like security, surveillance or human-computer interaction. Difficulties are due to appearance changes of the object and its background, motion variations or occlusion. Many methods have been proposed in the literature. We refer to the survey of Yilmaz [START_REF] Yilmaz | Object tracking: A survey[END_REF] for a comprehensive coverage of the problem. In this paper, we propose an original tracking method combining two algorithms: the Particle Filter and the Generalised Hough Transform. Isard adapted the particle filter for object tracking in his condensation algorithm [START_REF] Isard | Condensation: conditional density propagation for visual tracking[END_REF]. It belongs to the family of the Bayesian filters, and is based on the estimation of the posterior distribution of the object state using a set of weighted particles. Each particle corresponds to an hypothetical state of the object (typically describing its position and velocity), while the weight represents the importance of the particle. The particle filter is composed of three steps: (1) Observation, to measure how each particle fits the model, (2) Resample, to keep the best particles, and (3) Propagation, to diffuse the set of particles according to a dynamical model. In tracking problems, the observation step is usually based on object visual representations. Different representations have been used to model the target object: Isard [START_REF] Isard | Condensation: conditional density propagation for visual tracking[END_REF] models the borders of the object using parametric B-spline curves, Nummiaro [START_REF] Nummiaro | A color-based particle filter[END_REF] uses colour histograms. Pérez [START_REF] Pérez | Color-based probabilistic tracking[END_REF] splits the object in different parts, each one modelled by a colour histogram. Maggio [START_REF] Maggio | Combining Colour and Orientation for Adaptive Particle Filter-based Tracking[END_REF] uses smoothed histograms of gradient orientation for rotation robustness purposes, combined with colour histograms to form a particle filter using two features. Brasnett [START_REF] Brasnett | Sequential Monte Carlo tracking by fusing multiple cues in video sequences[END_REF] also uses multiple features (colour, edge and texture), and proposes a way to adaptively change their relative weights over time, to improve the performance of his tracker. The resample step is also an important process: Arulampalam compares different resampling methods in [START_REF] Arulampalam | A tutorial on particle filters for online nonlinear/non-Gaussian Bayesian tracking[END_REF]. Soto [START_REF] Arulampalam | Self adaptive particle filter[END_REF] proposes a method to adapt the number of particles according to the quality of the observation, and also a way to adapt the dynamical model used in the propagation step over time. The Hough Transform was originally proposed [START_REF] Duda | Use of the Hough transformation to detect lines and curves in pictures[END_REF] to detect lines and circles in binary images. Ballard [START_REF] Ballard | Generalizing the Hough transform to detect arbitrary shapes[END_REF] extended the algorithm to arbitrary objects, in the so-called Generalised Hough Transform (GHT). Originally, it consists in building a model using points sampled from the object boundary, and a function mapping the local orientation to the position of every point with respect to the model centre. Initially proposed as a detection algorithm, the GHT was extended to object categorisation [START_REF] Leibe | Combined object categorization and segmentation with an implicit shape model[END_REF]. In this last method, the representation was based on image patches built around interest points. The GHT was also used for object tracking. Duffner [START_REF] Duffner | PixelTrack: a fast adaptive algorithm for tracking non-rigid objects[END_REF] combined the GHT with a foreground/background model, whereas Godec [START_REF] Godec | Hough-based tracking of non-rigid objects[END_REF] enriched the model representation, using random forest theories. A similar model was proposed by Gall [START_REF] Gall | Class-specific hough forests for object detection[END_REF] for action recognition. The method proposed in this paper combines Hough-based tracking-by-detection with global and predictive particle filtering. The interest is to exploit the complementarity of the two approaches: the GHT is less sensitive than the particle filter to strong apparent accelerations. On the other hand, the GHT lacks any predictive capability and is more sensitive to object deformations. Moreover, while the particle filter is based on a global object model, the GHT uses a pixelbased representation. Finally, even in their simplest form, both algorithms have decent tracking performance, are very versatile, and can be improved in many ways. Our main contribution is the use of the output of the two algorithms to update the object model. It is based on the back-projection map from the GHT, which has proved a powerful tool. Razavi [START_REF] Razavi | Backprojection revisited: Scalable multi-view object detection and similarity metrics for detections[END_REF] showed that the back-projection can be a base for a multi-view GHT, and can measure the similarity between detections and training examples. It was also used to update the object model, by indicating if a pixel is likely to belong to the object. Our approach is different to those proposed by Duffner or Godec [START_REF] Duffner | PixelTrack: a fast adaptive algorithm for tracking non-rigid objects[END_REF], [START_REF] Godec | Hough-based tracking of non-rigid objects[END_REF], who use segmentation algorithms to strictly determine which pixels belong to the object. Our method uses two multivalued maps, concurrently produced by two approaches and combined to soflty evaluate the level of belonging to the object for each pixel. Combining a detector with a particle filter for object tracking has already proven effective: Breitenstein [START_REF] Breitenstein | Robust tracking-by-detection using a detector confidence particle filter[END_REF] proposed a multi-person tracker combining a detector (based on HOG [START_REF] Dalal | Histograms of oriented gradients for human detection[END_REF] or on GHT [START_REF] Leibe | Robust object detection with interleaved categorization and segmentation[END_REF]), a classifier and a particle filter. However, while the observation of his particle filter combines the detector and classifier outputs, in our case the observations are independent and only based on visual features. This paper is organised that way: Section II focuses on the particle filter, then Section III deals with the GHT. Section IV shows how the two algorithms are combined, and Section V presents results and discusses the performance of the tracker.

II. PARTICLE FILTER

A. Principle

The goal of the particle filter is to estimate the object state using a set of particles. In our case, we will work with a fixed number of particles: N p = 500. Each particle represents an hypothesis of the state, and is associated with a weight, whose role is to evaluate the adequation between the hypothesis and the model. In our experiments the position of the object is given by its bounding box. We then used the following fourdimensional state space:

• (x, y): coordinates of the bounding box centre • (w, h): width and height of the bounding box We will not consider the orientation of the object. The state of a particle i at frame t will be denoted by

X i t = {x i t , y i t , w i t , h i t }.
The three steps of the particle filter are defined as follows, for each frame:

• Observation: for each particle with state X i t , we associate a colour-based weight ω i t , obtained by comparing the colour histograms Ĥt of the object model and H i t of the image at frame t within the bounding box given by X i t . The RGB colour space is used, using n c × n c × n c bins (with n c = 8), and the histogram distance is based on the Bhattacharyya coefficient:

d c ( Ĥt , H i t ) = 1 - u∈{1,...,nc} 3 Ĥt (u) • H i t (u) (1)
and then, as proposed by Nummiaro [START_REF] Nummiaro | A color-based particle filter[END_REF], the colour weight is set as:

ω i t = exp [-λ • (d c ( Ĥt , H i t )) 2 ] (2)
with λ = 50.0, chosen experimentally. The use of colour histogram has the advantage to be robust to scaling and rotations issues.

From the observations, the global state of the object is estimated as:

Xt = Np i=1 X i t • ω i t ( 3
)
with N p the number of particles. • Resampling: Arumpamlam proposed [START_REF] Arulampalam | A tutorial on particle filters for online nonlinear/non-Gaussian Bayesian tracking[END_REF] different methods to resample the particles. For the sake of simplicity, the particle set was generated using a multivariate Gaussian Process. For each particle X i t = {x i t , y i t , w i t , h i t }, the multivariate Gaussian Process is done using four independent Gaussian Process :

-For the spatial states (x i t , y i t ), we use two independent Gaussian of parameters (x t-1 , c • ŵt-1 ) and (ŷ t-1 , c • ĥt-1 ) (x t-1 and ŷt-1 corresponding to the spatial estimation from the previous frame, while ŵt-1 and ĥt-1 corresponding to the width and height estimated in the previous frame). We choose c =

1 2•(2• √ 2•ln 2) , because (2 • √ 2 • ln 2)
is the coefficient linking the full width at half maximum (F W HM ) of a Gaussian and its standard deviation

σ: F W HM = (2 • √ 2 • ln 2) • σ.
To goal is to ensure that most of the particle are inside the last bounding box estimated.

-For the scale states, given the estimation from the last frame ( ŵt-1 , ĥt-1 ), we are generating a multiplicative coefficient β using a centralised Gaussian with a standard deviation of 0.05. Then

w i t = β • ŵt-1 and h i t = β • ĥt-1 .
• Propagation: by denoting A t a transition vector modelling a dynamical model, and W t a vector modelling a Gaussian noise, the propagation model to update the particle state from X i t+1 to X i t is given by:

X i t+1 = A i t • X i t + W i t (4)
The first matrix is used to describe the dynamical changes of the model (velocity, scale change), while the second is used to diffuse the particle set in the state space, but given our resampling model, the operation of diffusion has already been done before. In the experiments shown in this paper, the transition matrix A i t = I 4 is the identity, and W i t is a null-vector.

III. GENERALISED HOUGH TRANSFORM

A. Construction of the model

The GHT model is based on a so-called R-table T, built using a reference point R (the centre of mass of the object for example) and the gradient map. The R-table is indexed by orientation values, and each entry contains a list of couples ( u, ω), where u is a displacement, and ω its weight. Let M be the magnitude map, and O the quantified orientation of the gradient (with n o = 8 bins).

For all pixels P such that M(P ) is above a certain threshold (set to 10 in our experiments), let u P = -→ P R be the displacement of P relative to the centre, and Θ P the quantified orientation of the gradient at P , we add in T(Θ P ) the couple ( u P , M(P )). Like Duffner in [START_REF] Duffner | PixelTrack: a fast adaptive algorithm for tracking non-rigid objects[END_REF], for each table indexed by its orientation, we only store the D displacements with the highest weights (in our case D = 50). It has the advantage to limit the memory use and also significantly improve the speed (otherwise, each table can contain hundreds of displacements).

B. Detection

Given an image, the goal is to detect an object using the R-Table . The decision is made using a accumulation map (the Hough transform), denoted HT. The magnitude map M and the quantified orientation map O are computed. For each pixel P with M(P ) above the threshold, and orientation Θ P , we consider all the n Θ P entries {( u j , ω j )} j≤nΘ P of the R-Table T at the index Θ P . For each couple j, the accumulator HT(P + u j ) is incremented by the quantity ω j . The detector output for the best object position (if needed) may be given by P max = arg max P HT(P ).

C. Back-projection map

Under this form, the GHT has some weaknesses: in spite of a certain robustness to small deformations, it is very sensitive to scaling or rotating. Moreover, in a tracking context, it lacks a model updating process. Godec and Duffner proposed in [START_REF] Godec | Hough-based tracking of non-rigid objects[END_REF] and [START_REF] Duffner | PixelTrack: a fast adaptive algorithm for tracking non-rigid objects[END_REF] to use a back-projection map as a support for updating. In their case, after having detected the pixel P max with the highest value in the accumulation map, they build the back-projection map by determining which pixels voted for P max , for they are the more likely to belong to the object. Our approach uses a "softer" version of the back-projection map: each voting pixel is attributed the value of its best pixel in the accumulation map, more precisely:

BP(P ) = max j≤nΘ P HT(P + u j ) (5)
The map is then soft and denser compared to the one proposed by Duffner or Godec, however, a lot of pixels from the background may have a non-zero value. A threshold is not satisfying, since background regions similar to the object may have good votes in bad detection situation. Godec proposed to use the GrabCut segmentation algorithm from [START_REF] Rother | Grabcut: Interactive foreground extraction using iterated graph cuts[END_REF] to extract the foreground, while Duffner used a probability model to segment the foreground from the background. In our case, we do not used hard segmentation but combine the rough backprojection map with the output of the particle filter.

IV. FUSION

A. Combine the outputs

The goal is to exploit the output data from both the particle filter and the back-projection map to determine which pixels are more likely to belong to the object, and use them to update the two models. As written in Section III, the back-projection map provides local confidence measure related to the GHT model (the R-table). However, our "soft" back-projection map may also give good score to pixels far from the object. We then combine this value with the information provided by the particle filter. Let the observation map of the particle filter be the image OP defined from the position states of the particles, such that for all i ≤ N p , OP(x i t , y i t ) = ω i t , and OP = 0 elsewhere. Multiply the back-projection map BP with OP would be the natural combination. However, the latter is only a sparse estimation of the density map, and so, we need to smooth it. While fitting it with a joint probability distribution can give good results, it relies on a global estimation scheme that may be time-consuming. Instead we spread it using a morphological dilation, with an elliptic structuring element SE of size (

ŵ Np 1 3 , ĥ Np 1 3
), with ŵ and ĥ the estimated width and height of the object. Finally the fusion map B is defined as follows:

B = BP • δ SE (OP) (6)
B is then normalised to [0, 1] and all values inferior to = 0.05 are set to 0 (corresponding to discard pixels whose importance in the fusion map are less than 5% of the maximum value).

B. Update models

Given the normalised fusion map B, the two models are updated as follows:

• For the particle filter, we evaluate the colour histogram (denoted H B t ) from the sub-image defined by the bounding box of all pixels P such that B(P ) > 0. Then, we update the model according to this formula:

Ĥt+1 = α • Ĥt + (1 -α) • H B t (7)
with α = 0.9.

• to update the R-table, the magnitude map M and quantified orientation map O of the gradient are computed just like in the construction of the R-table. Then, let R be the centre of mass of B. For all non-zero pixels P from B, we consider its quantified orientation Θ P , and its fusion map value B(P ):

if u P = -→ P R is not in the table T(Θ P ), then we add the entry with a weight M (P ).

if it is already in, we update its weight ω Θ P , u P using both the magnitude M(P ) and the fusion map value B(P ):

ω Θ P , u P ← α • ω Θ P , u P + (1 -α) • B(P ) • M(P ) 2 (8)
To avoid overloading the memory and slowing down the process, the R-Table size is still limited to the D = 50 displacements with highest weights for each Θ P index.

C. Tracking decision

The final step is to provide the estimated state of the tracked object. In one hand, the particle filter gives and estimation of the object state. In the other hand, following the idea that a pixel in back-projection map is as high as it is more likely to belong to the object, we will use this back-projection map to estimate the state. Like the updating process, the use of the fusion map B is the natural combination. So, we decide to estimate the state of the tracked object by the bounding box of all non-zero pixels of B. The Figure 5 synthesises the complete tracking process as presented in this section and the two previous ones. 1. Generalised Hough Transform used alone. After several frames, the peak of the accumulation map (second image) becomes less visible. The third image represents the back-projection map, and the third map shows the thresholded version (using Otsu's method explained in [START_REF] Otsu | A threshold selection method from gray-level histograms[END_REF]). Some pixels from the background, in the bottom-right corner, are still visible Fig. 3. Success case for the combined tracker for position and scale changes. Alone, the particle filter was not able to detect the scale change. But with the back-projection map (the third image), pixels from the border of the mug, having high value, are still present in the fusion map (the fourth image), making success estimation of the scale and the position.

V. EXPERIMENTS

We have performed a series of preliminary experiments, whose purpose is to show that the combination of the two algorithms provide significant improvement with respect to each one used individually. In the working sequence, the objective is to track the mug. Basically, the movement is composed of a translation and a scaling, followed by a descaling. Both algorithms are used in their simplest possible form. On the one hand, the GHT is rapidly turning bad (See Figure 1): the accumulation map shows multiple local maxima with high values, and the lack of information from the previous frames (the GHT is just a detector) does not provide a way to select the best maximum. Moreover, without other information (such as the classification model used by Duffner [START_REF] Duffner | PixelTrack: a fast adaptive algorithm for tracking non-rigid objects[END_REF]), or an effective segmentation algorithm (like GrabCut used by Godec [START_REF] Godec | Hough-based tracking of non-rigid objects[END_REF]), extracting the object boundaries is a difficult task. On the other hand, while the Particle Filter, in the simple form proposed in our paper, gives decent results in the spatial estimation, the absence of propagation model makes it impossible to detect scale change (see Figure 2). However, it is possible to improve the tracking using differents solutions: using a propagative model taking in consideration the motion can be a solution. By combining the two algorithm, the tracker turns more effective to scaling issues (see Figure 3). However, it possesses the same weakness as the simplified Particle Filter: without a strong model of resampling and propagation, it is difficult to follow an abrupt variation of scale, and also possess the inherent problem of the GHT: without scale or orientation parameters, the detectors fails in rotation or scaling task (see Figure 4).

VI. CONCLUSION

We have presented an original and versatile object tracking algorithm combining Hough based tracking-by-detection and Particle Filtering. These two algorithms present fundamental complementarities: the GHT is insensitive to the acceleration while the Particle Filter has a strong predictive basis. Furthermore they operate on two different levels: pixel level for the GHT, and global level for the Particle Filter. The main rationale for this approach is in the hypothesis that the combination of two complementary methods under their most basic form provide more benefits that sophisticating each approach individually. Our contribution is focused on the method to fuse the outputs of the algorithms: combining the back-projection map from the GHT and the estimation states given by the Particle Filter provides an effective method to evaluate the level of belonging to the object at the pixel level. Our experiments showed that working with very light versions of the two algorithms, with a small R-Table for the GHT, and a fixed and small number of particles, the combined Fig. 4. Fail case for the combined tracker. In the first line, after several frame of fast descaling, the tracker is no longer able to estimate the scale. In the second line, after several frames of rotation, the GHT causes the fail: the second image shows different maximums at different places, and near of the center of the mug, several local maximums exists, perturbating the updating process, by adding datas from the background. Fig. 5. Summary of the tracker. At the top, the object is represented by its bounding box (the red rectangle) described using two features: gradient and colour. The goal is to estimate the state vector at each frame of the sequence. The GHT and the Particle Filter are operating independently: the first one generates the back-projection map, the second one the (spreaded) observation map. The two output are combined (Fusion step), while the tracker output is given by the Particle Filter.

tracker provides significant improvement compared to the two algorithms used independently. In the close future, we will focus on the extension of the state space to orientation parameter. We will also improve the fusion process, by proposing alternate combination between the two confidence maps. At middle-term, we will consider the different enhancement possibilities, at the visual representation level (e.g. use of different features) and on the algorithm itself (e.g. by implementing stronger resampling and propagation processes), with the objective to maintain the simplicity and lightness of the combined tracker. Quantitative evaluation of enhanced versions will be performed using academic datasets.

At long-term, the goal will be to extend the tracker to the multi-target context [START_REF] Breitenstein | Robust tracking-by-detection using a detector confidence particle filter[END_REF].

Fig.

  Fig.1. Generalised Hough Transform used alone. After several frames, the peak of the accumulation map (second image) becomes less visible. The third image represents the back-projection map, and the third map shows the thresholded version (using Otsu's method explained in[START_REF] Otsu | A threshold selection method from gray-level histograms[END_REF]). Some pixels from the background, in the bottom-right corner, are still visible

Fig. 2 .

 2 Fig. 2. Particle Filter used alone. Each coloured pixel corresponds to the projection of the particle in the spatial space. Its colour represents its weight (Equation2): red corresponding to high weight, blue to low weight. The rectangle is the state estimated with the Particle Filter, processed alone, using the equation 3. Its colour corresponds to the weight of the observation, for its state.

ACKNOWLEDGMENT

We gratefully acknowledge funding support by the French General Directorate for Armament (DGA).


Alaa Abi Haidar

From Inter-Annotation to Intra-Publication Inconsistency

What are effective ways to help people with chronic illness, e.g. diabetes and heart disease? Computational linguistics relies on human-annotated data to train machine learners. Inconsistency among the human annotators must be carefully managed (otherwise, the annotations are useless in computation). How can this annotation process be made scalable?

INTRODUCTION

Mathematical modeling and simulation help us understand the underlying mechanisms behind complex and barely understood systems, such as immune systems in order to advance biomedical and drug studies and cure diseases [10]. However, mathematical modeling requires huge amounts of parametric data, usually published in experimental and theoretical manuscripts and dispersed in the scientific literature. Pubmed comprises more than 19 million scientific articles [9] and this amount is growing at astounding rates. The manual extraction of valuable information and their classification into predefined labels, such as parametric values, units and species names, is very costly and inefficient. Hence, we use text mining, and more specifically, named entity recognition [8], in order to automatically and accurately extract and classify numerical and textual entities, that can later be plugged into mathematical models and simulated.

Nevertheless, a significant number of parametric values describing experimental results in the bio-literature are inconsistent due to variations in experimental approaches or imperfection in the experiments [19]. For example, the amount of T cells that mature from the thymus through the process of negative selection1 remains a huge debate [1][2][3][4][5][6][7] and may vary between "less than 5%" to "10%". Our study focuses on analysing all reported rates for such processes in biological and other complex systems in order to study their statistical variations while identifying average values and outliers.

BACKGROUND

Most techniques for named entity recognition rely on supervised machine learning that require human-annotated training data. Studies have shown that there is at least 25% human inter-annotation disagreement (inconsistency) when annotating biological named entities [24]. According to the G-theory [20], sources of inconsistency might be of external influences like alterations in the tools used for annotations, increasing time pressure, removal or adding of rewards, or changes in the annotation scheme. A study reports 55% and 82% F1 scores for exact and relaxed inter-annotation agreement respectively when the task was to extract 2 interactions between enzymes and marine drugs in over 230 full-text articles [START_REF] Rak | Database[END_REF]. In another study, the relaxed inter-annotator agreement showed that 94% of the time curators were precisely extracting GO annotation from the literature and 3 72% of the time curators recalled all possible valid GO terms from the text [14]. Yet another study reports an inter-annotator agreement rate of over 60% for triggers and of over 80% for arguments using an exact match constraint [16]. More recently, a study reported inter-annotator agreement (IAA) F-measures for medication names and medication types, 94.2% and 88.2% respectively [START_REF] Zhai | [END_REF]. However, the rates of inter-annotation agreement vary not only from study to study but also from domain to domain. For instance, Wiebe et al. report 82.0 F1 score of human annotation agreement for opinion expression [START_REF] Wiebe | Language Resources and Evaluation[END_REF].

Consequently, the human inter-annotation inconsistency creates a gray area of uncertainty that the machine learner depends on to create fine-tuned rules and exceptions. Furthermore, human annotation is often used as the gold standard for evaluating machine learning methods [15] and therefore it is very important to have as few disagreements as possible. Nevertheless, the human inter-annotation disagreement can be reduced by using strictly agreed-upon annotations, the reasonings of a single annotator, majority rules or by identifying mislabeled annotations [17,18].

RESULTS AND DISCUSSION

Here, we attempt to quantify human annotation inconsistency based on a biological article annotated by several annotators. More specifically, we study the annotation inconsistency of a biomedical article [24] that is annotated by three experts for 9 categories of named entities. 1. The number of annotated entities in a biomedical article [24] for each of the 9 categories by each of the three annotators

The categories describe biological concepts. In the following example "The T cell proliferation is 0.4 cells/hr" can be annotated as follows: "T cell proliferation" describes a PROCESS, "0.4" a NUM (numerical value) and "cells/hr" a UNIT. More complex entities are harder to classify into concepts which may create annotation inconsistencies between several annotators.

For the example at hand, we identify inconsistencies varying from 1.5% to 8.3% out of a total of 10052 terms as shown in table 2. Our average values are below those reported by [24]. However, that might be due to the fact of working on a different dataset and with different entities. Next, we study to what extent inter-annotation inconsistency in training data can influence the robustness of machine learning and the predicted results. We hypothesize that factors such as the size of annotated training data, the number of annotators, the number of class labels, and the over-fitness of supervised machine learners play major roles in the robustness of the learning and the classification results. We attempt to answer these questions using empirical approaches inspired by [21] a study of robustness when classifying noisy land rover data from satellite images.

Al

EVALUATING INTER ANNOTATOR AGREEMENT

Determining inter annotator agreement (IAA) can be a daunting task for reasons such as proper choice of agreement coefficients and lack of consensus on the interpretation of such coefficients [26].

Taking into account the recommendations of Artstein and Poesio [26 -590] considering the better quality measure given by chance corrected coefficients as compared to simple percentage agreement, we have performed the reliability testing over the annotations performed over our corpus of text. Given that we are in a scenario of multi-annotators (namely, three) using a nominal variable, we will discuss the results we've obtained for the computation of the adapted versions of Cohen's ᵰ and Fleiss' multi-ᵰ as well as for Krippendorff's ᵯ.

As a reminder, the first 2 coefficients above are based on the basic coefficients used in 2-annotators scenarios: Scott's ᵰ (1955), Cohen's ᵰ (1960).

,

π, k = 1 -A e A -A o e
where A o is the observed agreement and A e is the expected agreement. [26 -559] As explained in [26 -560], the difference between ᵰ and ᵰ lies in the assumptions made to compute the probability for a coder (annotator) to categorize an utterance in a certain category.

The indices we used are the generalizations of Scott's ᵰ and Cohen's ᵰ made by Fleiss (1971) and Davies and Fleiss (1982) respectively.

To generalize, Fleiss' multi-ᵰ uses a different interpretation of the observed annotation A o , namely the pairwise agreement, which is the number of pairs agreeing on an utterance out of the total number of pairs of coders. Equally, multi-ᵰ involves the computation of the expected agreement A e based on individual coder marginals.

Finally, Krippendorff's ᵯ is a versatile coefficient which addresses the limitations of (multi-)ᵰ and (multi-)ᵰ regarding the equal treatment of all disagreements [26 -564]

1 α = -D e D o
, where D o is the observed disagreement and D e is the expected disagreement. [26 -565, 566] We computed the coefficients in the following 2 situations [START_REF] Freelon | ReCal3: Reliability for 3+ Coders[END_REF][START_REF] Freelon | [END_REF] Looking at the "classical" agreement reporting, the average pairwise percentage agreement is at 94.2%, with the best agreement rate being between coder 1 and coder 2 at 98.5%.

As expected, the values for multi-ᵰ and average pairwise multi-ᵰ are approximately equal and since all disagreements are treated equally, Krippendorff's ᵯ is also nearly equal to the two before. The approximate value of 0.82 for multi-ᵰ, classifies the annotation process as "perfect" according to the strength scale given by Landis and Koch (1977) [26 -576]. In this case, almost 8000 utterances which have been all annotated as Other by all three coders were completely discarded from the reliability study.

The classic pairwise percentage agreement thus records a drop of the average value (75%); however, it is interesting to notice that the agreement between coders 1 and 2 stays as high as in the original situation (93.4%) which indicates a very high agreement for their annotations on the "main" classes (excepting the all-O(ther)) as well.

As before, the values of the three coefficients are (not surprisingly) almost equal, but are dropping this time, somewhere around 0.67, which on the strength scale of Landis and Koch (1977) [26 -576] indicates only a "substantial" agreement in the annotation process, making the annotation suitable for "tentative conclusions". Interesting enough, this is exactly the value which was set as the original threshold by Krippendorff (although he referred to it as "highly tentative and cautious") before he reviewed it and later set it at 0.8.

CONCLUSION

In this manuscript, we identify and discuss two forms of inconsistencies, one in published parametric results that we are studying in collaboration with immunologists , and another in manually human annotated data for machine 4 learning. Both forms of inconsistencies influence the accuracy of biomedical research to a significant extent that we are interested in quantifying in our research. We expect our study to shed a light on both forms of inconsistencies, ones resulting from human inter-annotation and those published in biological literature.

ACKNOWLEDGEMENT

"This work has been done within the LABEX OBVIL project, and received financial state aid managed by the Agence Nationale de la Recherche, as part of the programme "Investissements d'avenir" under the reference ANR-11-IDEX-0004-02". We thank the expertise and effort of Veronique THOMAS-VASLIN and Florence GIESEN for the help with the annotation and valuable discussions.

BIBLIOGRAPHY

: 1 .

 1 Full annotations, taking into account the large number of utterances classified as O(ther) by all annotators (Coders = 3, Utterances = 10043).

2 .

 2 Altered annotations, discarding the utterances classified as O(ther) by all annotators (Coders = 3, Utterances = 2349).

Table

  The following table lists the number of annotated entities for each of the 9 categories by each of the three annotators:

	almeida_annotation_	almeida_annotation_Veroni	almeida_annotation_Floren
	Al.tag	que.tag	ce.tag
	17 UNIT	17 UNIT	28 LOCATION
	22 LOCATION	20 LOCATION	30 UNIT
	34 NUM	65 NUM	38 INDIVIDUAL
	63 INDIVIDUAL	74 INDIVIDUAL	58 NUM
	227 METHOD	226 METHOD	161 METHOD
	255 COFACTOR	231 COFACTOR	356 COFACTOR
	344 PROCESS	360 PROCESS	454 PROCESS
	950 POPULATION	907 POPULATION	904 POPULATION
	8140 O	8152 O	8014 O

Table 2 .

 2 Inter-annotator inconsistencies varying from 1.5% to 8.3% out of a total of 10052 terms

	Ver	Flo

Table 3a

 3a

	multi-ᵰ	Ao		Ae
	0.827	0.942	0.662
	Table 3c. multi-ᵰ (based on Scott's ᵰ)
	Krippendorff's ᵯ	No. of decisions
	0.827	30129
	Table 3d. Krippendorff's ᵯ
	. Average Pairwise Agreement [%]
	Avg. Pairwise CK	Pairwise (1-2) CK	Pairwise (1-3) CK	Pairwise (2-3) CK
	0.828	0.953	0.774	0.758

Table 3b. Average Pairwise Multi-ᵰ (based on Cohen's ᵰ)

Table 4a

 4a

	Krippendorff's ᵯ	No. of decisions
	0.674	7047
	Table 4d. Krippendorff's ᵯ
	. Average Pairwise Agreement [%]
	Avg. Pairwise CK	Pairwise (1-2) CK	Pairwise (1-3) CK	Pairwise (2-3) CK
	0.675	0.914	0.569	0.543
	Table 4b. Average Pairwise Multi-ᵰ (based on Cohen's ᵰ)
	multi-ᵰ	Ao		Ae
	0.674	0.75	0.233
	Table 4c. multi-ᵰ (based on Scott's ᵰ)

T cells that recognize self antigens are eliminated in the thymus through a process known as negative selection so that they do not bind to self and cause autoimmune diseases

Unlike exact matches, relaxed matches may span over less or more words to describe the same concept.

The Gene Ontology, or GO, is a bioinformatics attempt to unify the representation of genes and gene product attributes across various species

The team of Integrative Immunology (I2) is based at the Pitié Salpêtrière. 83, Bld de l'hôpital. 75013 Paris, France