Irving O Morales

E Landa

P Stránský

A Frank

An improved unfolding by detrending of statistical fluctuations in quantum spectra

Keywords: numbers: 05.40.-a, 05.45.-a, 05.45.Tp

A fundamental relation exists between the statistical properties of the fluctuations of the energy level spectrum of a Hamiltonian and the chaotic properties of the physical system it describes. This relationship has been addressed previously as a signature of chaos in quantum dynamical systems. In order to properly analyze these fluctuations, however, it is necessary to separate them from the general tendency, namely, its secular part. Unfortunately this process, called unfolding, is not trivial and can lead to erroneous conclusions about the chaoticity of a system. In this paper we propose a technique to improve the unfolding procedure for the purpose of minimizing the dependence on the particular procedure. This technique is based on detrending the fluctuations of the unfolded spectra through the empirical mode decomposition method.

I. INTRODUCTION

The relation between classical dynamical systems and their quantum counterparts is an important aspect of contemporary research. Of particular relevance is the study of the so called quantum chaos, the relation between classical chaos and quantum mechanical systems. One of the fundamental signatures of quantum chaos is the conjectured link between the statistical fluctuations of the energy spectrum and the integrability or chaotic properties of the Hamiltonian [1][2][START_REF] Stöckmann | Quantum Chaos[END_REF]. This relation has been studied in a large diversity of systems, including microwave cavities, atoms, molecules and nuclei. This kind of spectral analysis has become one of the main tools for the study of quantum chaos, whose fundamental objective is to characterize the properties of these fluctuations [START_REF] Brody | [END_REF] To tackle this problem it is assumed that the density of states ρ(E) can be separated into a secular and a fluctuating part

ρ(E) = ρ(E) + ρ(E). (1)
While the secular part ρ(E) usually indicates particular characteristics of the system under consideration, the fluctuating part reflects universal properties of physical systems. However, to separate these two components of the spectral density is a major problem. The procedure used to extract the information of the fluctuating part in the spectra is called unfolding. In essence, this procedure consists in mapping the spectrum of the excited states E i into a dimensionless spectrum ǫ i , whose main characteristic is a mean level density normalized to 1. It has been widely appreciated that this process is not trivial [START_REF] Brody | [END_REF], but the unfolding details are usually omitted in the literature. Of special importance to the present work is the analysis reported on [5], which demonstrates that the statistics that measures long-range correlations strongly depends on the unfolding procedure and in the end, this dependence can give rise to misleading results in regard to the chaoticity of quantum systems. The main goal of the present paper is to propose an unfolding technique which leads to reliable signatures of chaos. This is achieved subtracting the residual trend which remains after the initial unfolding process. This is achieved by means of the empirical mode decomposition technique.

II. UNFOLDING PROCEDURE

Given a spectral density ρ it is possible to define an accumulated level density N (E), which measures the number of levels up to a certain energy E in the system

N (E) = E -∞ dE ′ ρ(E ′ ). (2)
In this way we can distinguish between a smooth and a fluctuating part

N (E) = N (E) + N (E). (3)
The unfolding procedures usually involves an estimation of a functional shape for the smooth part of the spectral density ρ(E) or for the smooth part of the accumulated level density N (E). Once the estimated accumulated density is known a new dimensionless spectrum ǫ i is obtained, where the levels are given by

ǫ i = N (E i ). ( 4
)
The problem is how to find an appropriate smooth estimation for the secular part of the level density. When possible, the spectral density is obtained theoretically using the knowledge of the low order energy moments, which are related to the trace of powers of the hamiltonian matrix [6]. If this is not possible, then it is necessary to obtain it from the spectrum itself. Several procedures have been proposed in order to unfold this kind of spectra. The easiest unfolding approach is to assume a polynomial dependence for the accumulated level density [7]. In this approach a high order polynomial interpolation for N (E) is taken, the only parameter being the order of the polynomial approximation. However, this estimation can be quite inaccurate, since the polynomial function may have no relation with the physical properties of the system. It is also possible to consider a linear fit of N (E) on an interval containing ν levels around the level E i , and a moving average then taken through the whole spectrum. The ν value is used as a free parameter [8]. More sophisticated approaches for the unfolding involve Fourier broadening of the step functions, which reflect the accumulated density or a Gaussian broadening of the level density. While usually this density is conformed by a summation of delta functions for every level, under this approach the delta functions are substituted by an average gaussian level density. It is also possible to include local effects in this procedure [9].

Regardless of which method is used to unfold the fluctuating part of the spectra and even if a good estimation for the accumulated level density is obtained, this approximation may not be good enough. If this is the case the unfolding does not capture the full behavior of the smooth part of the spectrum. The net effect is that spurious long-range correlations are introduced. These spurious correlations are responsible for the misleading signatures of quantum chaos that may occur [5]. The technique that we propose here does not depend on apriori selection of unfolding functions and, as we shall show, significantly improves the estimation for the smooth part. Throughout this work we carry out an initial unfolding with polynomial approximations as an example, but the results do not depend on this choice and it is possible to improve the estimation starting with any of the previously mentioned procedures.

Once a dimensionless spectrum is obtained, we can calculate the nearest neighbor level spacings defined by

s i = ǫ i+1 -ǫ i , i = 1, ..., N -1. (5)
Recently, a proposal was made to study the fluctuating part of a quantum spectrum as a time series, which gives rise to an alternative signature of quantum chaos [10]. This approach characterizes the spectral fluctuations with the δ n statistic [START_REF] Mehta | Random Matrices[END_REF] defined by

δ n = n i=1 (s i -s ). ( 6
)
The discrete function δ n measures the deviations between the energy of the (n + 1)-th unfolded state and the corresponding energy of an equally spaced energy sequence. The main objective of the analysis is to propose a correspondence between the δ n statistic and time series, by equating time to the index n, which represents the ordering number of the energy levels. It is then possible to apply the usual techniques of time series analysis. As we mentioned above, the result of a defective unfolding is the introduction of spurious correlations in the spectrum fluctuations. In a properly unfolded spectrum all traces of the secular component of the signal should in principle be removed. In terms of time series, this correlations are represented as global trends on the δ n statistic, i.e. a smoothly varying behavior is present. If the energy spectrum is correctly unfolded then the δ n should consist only of fluctuations around its constant mean value. Figure 1 shows the effect of unfolding with polynomials of different order on the same energy spectrum. The spectrum was obtained from a gaussian orthogonal ensemble (GOE) matrix, a system which is regarded as a paradigm for the study of quantum chaos [START_REF] Stöckmann | Quantum Chaos[END_REF]. The energy spectrum was unfolded with three different polynomials and the resulting δ n for each case is shown in figure 1. The bottom (red online) signal shows a δ n obtained with a fifth order fit, the middle (blue online) one shows the one obtained with a ninth order fit and the top (green online) signal shows the δ n obtained for a fit with a polynomial of order 13. The plots were displaced along the vertical axis to avoid overlapping. The dashed lines show the mean for each case. If the unfolding procedure were correct, the time series should fluctuate around this mean, but a clear global trend is still present. This trend is shown as a black solid line on each case. In order to solve the problem of an incorrect unfolding we propose to detrend the δ n time series using the EMD technique.

III. TRENDS AND EMPIRICAL MODE DECOMPOSITION

The trend is a basic concept in a large amount of subjects in physics, economics, statistics, sociology, etc. The concept is used specifically in the analysis of data sets. While an intuitive idea of the trend can be used, there is no precise definition of it. In most of the data analysis applications, the trend dominates the global behavior of the data. If a more precise analysis of the phenomena involved is required, then the trend must be substracted in order to reveal the relevant details. Unfortunately, a general algorithm to detrend a time series does not exist, especially in the case where the time series is nonlinear and nonstationary. As in the case of the unfolding procedure, an accurate and detailed description of the detrending procedure is usually avoided in the literature. Recently an attempt to mathematically define the trend has been proposed [START_REF] Wu | [END_REF], where the authors put special emphasis on the fact that the trend is part of the data, an intrinsic property of it. Therefore, in order to extract the trend an adaptive method is required, namely, the trend must be derived directly from the data. The empirical mode decomposition technique [13] is especially suited to accomplish a detrending process fully based on the data.

The empirical mode decomposition (EMD) proposes the expansion and decomposition of the data in terms of intrinsic mode functions (IMF) which are derived from the data and conforms an adaptive local and nonlinear basis. This basis is constructed using the physical time scales that characterize the oscillations of the phenomena. The EMD assumes that the time series under analysis is composed of the superposition of several nonstationary processes that take place in different time scales.

The essence of the method is to identify the oscillation modes of the data in a characteristic time scale and then decompose the data into these modes of oscillation. The local time scale of an oscillation is identified by the time lapse between the successive alternations of local maxima and minima. The time series is composed of several oscillations mounted on oscillations corresponding to longer time scales, which in turn are also riding bigger oscillations, and so on. The process of sifting, described below, has been designed in order to untangle these oscillations and extract them from the time series. The oscillating functions in which the time series is decomposed are called intrinsic mode functions. These intrinsic mode functions are defined as functions in which a) the number of extrema and the number of zero crossings are the same or differ at most by one, and b) at any point the mean value of the envelope defined by the local maxima and the envelope defined by the local minima is zero. In general the data oscillations involve many time scales and it is necessary to decompose the time series in its IMFs. The process of sifting to decompose the data into IMFs is as follows:

1. Localize all the local maxima in the time series and connect them with a cubic spline, which defines the upper envelope. Repeat the procedure with the local minima defining the lower envelope.

2. Calculate the mean value of the upper and lower envelopes m (1) , as well as the difference with the data points

x j -m (1)
j = h (1) j , (7)
Since the envelopes include all the extrema, the oscillation between them, h (1) , is the oscillation of the data with the finest time scale. Ideally this oscillation should be an IMF, but in practice the process produces new extrema and more iterations are needed to extract an IMF. In this new iterations the obtained component is used as the new data set.

3. Repeat steps 1 and 2 using the newly obtained component as data

h (1) j -m (1,1) j = h (1,1) j . ( 8
)
Stop when the last component is an IMF, i.e., when all local maxima are positive and all local minima are negative, and the mean of the envelopes is locally zero, which ensures that the wave is symmetric. The IMF is obtained after k iterations

h (1,k-1) j -m (1,k) j = h (1,k) j . (9)
The function h (1,k) is then designated as an intrinsic mode function of the data and as a component of our new basis

c (1) = h (1,k) . ( 10
)
4. Subtract the obtained IMF from the original data

x j -c (1)
j = r (1) j . (11)
The residue r (1) is a new version of the data where the finest oscillations have been removed. However, it still contains information of longer period components, so it is necessary to remove these additional components.

5. Repeat steps 1 to 4, in order to extract all the IMF components in the original data, using as new data on each iteration the residue obtained in the previous one

r (1) -c (2) = r (2) , ..., r (m-1) -c (m) = r (m) , (12)
where m is the number of IMFs in the data. This process can be stopped when the mth residue is a monotonic function and no more IMFs can be extracted. The last residue is called the trend of the data. It is important to notice that any residue constitutes a trend for the previously extracted oscillation, i.e., r (i) is the trend followed by the c (i) intrinsic mode function. For a discussion about the specific details of the algorithm and details about its convergence see [13][START_REF] Huang | Hilbert-Huang Transform and its applications[END_REF][15].

After following this procedure, it is possible to express the original data in terms of the IMFs as

x j = m i=1 c (i) + r (m) . ( 13
)
We are interested in the last residue r (m) , which is the global trend.

The trend should be a constant for a δ n obtained from a correctly unfolded spectrum. Figure 1 shows as solid black lines the trends obtained through the EMD for the same spectrum, previously unfolded with polynomials of order 5, 9 and 13. Two undesired effects are apparent in this plot, the first one related to the global trend, which is still visible after the unfolding procedure. In each case it is evident that the EMD technique captures the global behavior, so it is possible to subtract it. The second effect is more subtle and has to do with the fact that the values for the δ n on the edges of the spectrum change abruptly, especially when the order of the original polynomial fit is low. If this is the case, the EMD decomposition is not well suited for these areas, due to the border conditions of the cubic spline interpolations. In order to avoid this effect we drop 5% of the eigenvalues at the beginning and at the end of the spectrum.

In order to decrease the dependence on the unfolding procedure, we propose to construct a new detrended fluctuation, subtracting the trend obtained with the EMD, δ ′ n = δ n -r (m) . Figure 2 shows the result of excluding 5% of the eigenvalues on the edges of the spectrum and subtracting the residual trend from the original δ n statistic (Fig. 1). It is remarkable how in a qualitative examination the main features of the fluctuations are present on each plot and apparently display the same features. A subtle difference is appreciated at the last part of the δ ′ n obtained with a fifth order polynomial. This plot has an extra oscillation around n = 1200 which is present due to the effect mentioned before, the abrupt change of δ n values on the edges of the spectrum. In this case dropping 5% of the eigenvalues may not be enough. However this case is useful to appreciate how the affected part is restricted to the end of the δ ′ n statistic, i.e., it is a local phenomenon which can be removed. To quantitatively show that the new δ ′ n statistics are independent of the unfolding is the main goal of the next section.

IV. MEASUREMENTS OF QUANTUM CHAOS

In order to demonstrate that the δ ′ n statistic, obtained subtracting the EMD trend, leads to chaos measurements independent of the unfolding procedure, we performed a numerical experiment. We constructed and diagonalized 100 GOE matrices with 2000 eigenvalues per matrix. Then the δ n and δ ′ n statistics were computed for each spectrum. Given these two statistics the ensemble behavior of the usual measurements of quantum chaos is compared for different orders of the polynomial fit. The GOE system was chosen because it is a well known example of a system presenting quantum chaos and, in addition, the theoretical result for the usual measurements of chaos are known [START_REF] Stöckmann | Quantum Chaos[END_REF]10,[START_REF] Mehta | Random Matrices[END_REF].

As a first test we choose a quantum chaos measure which is computed using the δ n statistic, the power law behavior of the power spectrum of the fluctuation time series [10]. The power spectrum of a given time series, x n , is defined as

S(k) = | Xk | 2 = 1 √ N n x n exp -2πikn N , ( 14
)
where Xk is the Fourier transform of x n and N is the

′ n = δ n -r (m)
where r (m) is the trend obtained through EMD. The spectrum was unfolded with polynomials of order 5 (bottom, red online), 9 (middle, blue online) and 13 (top, green online). The dashed lines show the mean on each case. Each plot was moved along the vertical axis to avoid overlapping.

the number of elements in the series. It has been demonstrated that the power spectrum of the fluctuations of quantum spectra follows a power law dependence S(k) ∼ 1/k β , this relation is only an approximation, valid without taking into account the highest frequencies of the spectrum. Therefore in the numerical analysis, we take only the region of frequencies which follows a power law behavior, [10,[16][17][18]. In the case of an integrable system the value of the exponent approaches β = 2, because the inter-level spacings behave as random independent variables [2]. For the GOE system the value β = 1 corresponds to the fact that correlations exists on all available scales. According to this, the power spectrum of a correctly unfolded GOE spectrum must follow the power law S(k) ∼ 1/k. The numerical experiment compares the effectiveness of detecting the presence of quantum chaos using both statistics, the δ n fluctuation obtained by the traditional unfolding of the spectrum and the one without trend, the δ ′ n . In order to compare the ensemble behavior of the system, an ensemble average power spectrum is computed using the 100 individual power spectra, which lowers dispersion and brings out the power law behavior of the system. The results of the numerical experiment are compared in figure 3. This figure shows the values of the exponent β obtained for different orders of the unfolding polynomial fit. The filled squares (red online) corresponds to the case of the average power spectrum of the δ n fluctuation and the empty circles (blue online) shows the case of the power spectrum of the detrended fluctuation δ ′ n . The error bars show the standard deviation using individual fits to the power spectra of each member of the ensemble. The gray dashed line indicates the theoretical value of the β exponent for the GOE system. It is evident that the value of β depends on the order of the polynomial used to unfold the quantum spectra. If the order of the polynomial is too low then spurious longrange correlations are introduced and the low frequency part of the power spectrum is enhanced producing power law exponents larger than β = 1. However, the use of the EMD to remove the residual trend eliminates these spurious correlations and improves the value of the exponent which gets closer to the theoretical value. When the trend is removed, for all the degrees of the polynomial fit lower than 15, the β value agrees with the theoretical value, except for the example corresponding to the third order polynomial, which is an extreme case.

As the order of the polynomial increases the β value approaches the theoretical one. If the order of the polynomial is big enough then the polynomial unfolding destroys correlations and the β value tends towards values lower than β = 1. This situation cannot be improved by the use of the EMD because important correlations are missing and cannot be recovered.

Using the detrended fluctuation δ ′

n it is possible to ob-tain a detrended dimensionless spectrum ǫ ′ i as

ǫ ′ i+1 = δ ′ i + ǫ 1 + i s , (15)
where δ ′ n are the fluctuations or the difference between the dimensionless quantum spectrum of the system ǫ ′ and a regular version of it, ǫ 1 + n s . Given this detrended spectrum it is possible to explore the dependence of other measurements of chaos as a function of the order of the polynomial used in the unfolding. The short-range correlations are usually analyzed by means of the nearest neighbor spacing distribution P (s). It is known that the spectra from GOE follow the Wigner surmise. This distribution is not sensitive to the unfolding details [5] and we only mention here that the Wigner distribution is conserved after removing the trend of the fluctuations.

The study of the long-range correlations in quantum spectra is traditionally done by means of the rigidity of the spectrum, which is represented by the ∆ 3 statistic. The ∆ 3 (L) is defined as the average over γ of the function:

∆ 3 (L; γ) = 1 L min a,b γ+L γ [N (E) -aE -b] dE, (16)
where N (E) is the integrated level density of the unfolded spectrum and a and b are parameters of a linear model which fits N (E) in the energy range (γ, γ + L).

Theoretically the values of the ∆ 3 (L) for the integrable system and for the GOE are well known [19], in the case of GOE ∆ 3 (L) = (Log(L) -0.0687)/π 2 while in the case of an integrable system ∆ 3 (L) = 1/15.

It is possible to compute the ∆ 3 statistic for the detrended spectrum and compare it with the ∆ 3 obtained from the traditionally-unfolded spectrum. Figure 4 shows this comparison. The upper part shows the ensemble average of the ∆ 3 obtained for the unfolded spectrum using different degrees for the polynomial fit. The dashed gray line shows the predicted theoretical behavior for GOE. We observe that for all polynomial degrees the expected result is followed up to L = 200, except for the third order polynomial which shows some small deviations. For greater values of L the low order polynomial spectra deviate systematically from the theoretical behavior due to the spurious correlations, which correspond to the residual trend. Since the trend can be represented as low frequency oscillations it is expected that the effect will be reflected in the long-range correlations region. The high order polynomial spectrum also deviates because of the effect mentioned before for the case of the Fourier spectrum analysis, i.e., the polynomial unfolding destroys important correlations. The lower part shows the ensemble average of the ∆ 3 computed using the detrended spectra from eq (15). The effect of subtracting the residual trend of the fluctuations is evident. In the region of small L the deviations for the third order polynomial are gone while in the region of larger L the low order polynomial behavior improves and follows the expected theoretical behavior. The case for the high order polynomials does not improve, because once important correlations are destroyed they cannot be recovered.

The ∆ 3 statistic is a good measure of the long-range correlations of a quantum spectrum, but it is difficult to obtain a quantitative deviation from the theoretical behavior, as done in the previous example of the Fourier power spectrum by means of the value of the exponent β. Recently another measurement of the chaoticity of a quantum system has been proposed [START_REF] Santhanam | [END_REF]. This method uses a detrended fluctuation analysis technique [21] (DFA) which has been directly related to the ∆ 3 [START_REF] Santhanam | [END_REF]. The DFA technique was proposed in the context of time series analysis as a method to investigate the long-range correlation properties of nonlinear time series. In order to apply the DFA the time series is integrated and then divided into nonoverlaping N/l boxes of equal length containing l data points. Then a local trend is defined for each box, i.e. a linear fit to the integrated time series. Then the variance between the integrated time series and the local fit is computed for each box and averaged over all the boxes of size l. The average variance F (l) depends on the box size. The analysis of this dependence allows to study the scaling properties of the system and the presence of long-range correlations. A linear relationship between Log(F (l)) and Log(l) will indicate a scaling behavior. The slope α in the Log(F (l)) vs Log(l) plot characterizes the scaling properties of the time series because it reflects a power law of the form F (l) ∼ l α . In the case of a regular uncorrelated system the theoretical value of the slope is α = 1/2. In the case of GOE the expected value is α = 1 [START_REF] Santhanam | [END_REF]. The DFA technique allows us to compare the δ n fluctuation with the detrended one, δ ′ n . A qualitative comparison on the DFA relation between Log(F (l)) and Log(l) for the δ n and δ ′ n fluctuations is shown on figure 5. Once again, the improvement due to the detrended process is noticeable. The upper part shows Log(F (l)) as a function of Log(l) obtained from a DFA calculation on the unfolded fluctuations δ n for different degrees of the polynomial. The theoretical value for GOE (α = 1) is shown as a dashed gray line. It is evident how for low order polynomials the result diverges from the expected theoretical behavior when l increases. Large values of l correspond with long-range correlations on the time series. This plot shows how if the degree of the polynomial is too low the unfolding is capable of representing the correlations correctly, up to a certain range. On the other hand, if the polynomial degree is too big then the unfolding procedure destroys important long-range correlations. This effect is also observable on the upper part of figure 5 but is more subtle. The lower part of the same figure shows the result of the DFA on the detrended δ ′ n fluctuation. It is clear how the divergence observed for the low order polynomials has decreased. All the cases have a linear behavior which corresponds to the theoretical case. For very large values of l the third order polynomial case begins to diverge, which means that some very long-range correlations are still missing. However, the range for which the DFA variance Log(F (l)) follows a linear relation has been substantially improved. The results shown in figure 5 demonstrate qualitatively how detrending the fluctuation improves the range for which the long-range correlations of the system can be detected. However, the DFA allows us to observe this improvement in a quantitative way through the calculation of the value of the α exponent. Figure 6 compares the values obtained for α as a function of the order of the unfolding polynomial using the fluctuation δ n (filled squares, red online) and the values obtained using the detrended fluctuation δ ′ n (empty circles, blue online). The theoretical behavior is shown as a gray dashed line. The improvement for small values of the polynomial degree is evident, after detrending the α exponent coincides with the predicted one. This improvement is observed up to the ninth order polynomial. For cases with larger values of the polynomial degree there is no improvement, again due to the effect mentioned before, i.e., a high order polynomial unfolding destroys important correlations which cannot be recovered detrending the fluctuations.

V. CONCLUSIONS

In order to study the statistical fluctuations of quantum spectra it is necessary to split the spectrum into two parts, a secular and a fluctuating one. This separation is done by means of an unfolding procedure. Several methodologies exist to unfold quantum spectra, but this procedure is not trivial. In this work we demonstrate that the usual measurements of quantum chaos can strongly depend on the unfolding procedure, i.e., an incorrect separation of the fluctuating part from the regular one can lead to erroneous conclusions about the chaoticity of a quantum system. In order to decrease this dependence we propose to remove the global trend which remains after an incorrect unfolding procedure. The detection of this global trend is the essential part of our proposal, which makes use of the Empirical Mode Decomposition techinique. Finally, we show that the detrending operation on the fluctuations removes to a large extent the dependence on the unfolding procedure for the usual measurements of the presence of quantum chaos. The effect of an incorrect unfolding is either to introduce spurious long-range correlations or to destroy important ones. In the first case the use of the EMD to detrend the fluctuations improves the range on which the correlations are correctly represented.

FIG. 1 :

 1 FIG.1: Color online. δ n statistic for a GOE spectrum unfolded with polynomials of order 5 (bottom, red online), 9 (middle, blue online) and 13 (top, green online). The dashed lines show the mean on each case and the black solid line shows the trend obtained trough the EMD. Each plot was moved along the vertical axis to avoid overlapping.

FIG. 2 :

 2 FIG.2: Color online. Detrended δ n statistic for a GOE spectrum, definied as δ ′ n = δ n -r (m) where r (m) is the trend obtained through EMD. The spectrum was unfolded with polynomials of order 5 (bottom, red online), 9 (middle, blue online) and 13 (top, green online). The dashed lines show the mean on each case. Each plot was moved along the vertical axis to avoid overlapping.

FIG. 3 :

 3 FIG. 3: Color online. Value of the β exponent of the power law behavior in the average power spectrum of the fluctuations as a function of the degree of the polynomial fit used to unfold the quantum spectra. Two cases are shown, the β exponent values for the δ n fluctuation (filled squares, red online) and the case for the δ ′ n fluctuation obtained detrending the previous one (empty circles, blue online). The dashed line shows the theoretical result for the GOE. The error bars show the standard deviation of the ensemble.

FIG. 4 :

 4 FIG. 4: Color online. Value of the ∆ 3 (L) averaged over the ensemble. The upper part shows the ∆ 3 computed using polynomials of different order to unfold the spectra, where the plots have been displaced to avoid overlap. The lower part shows the ∆ 3 calculated from detrended spectra obtained from the detrended fluctuation δ ′ n . Plots for different degrees of the polynomial fit used to unfold the quantum spectra are shown in both cases. The dashed line shows the theoretical result for the GOE.

FIG. 5 :

 5 FIG. 5: Color online. Value of DFA variance Log(F (l)) averaged over the ensemble. The upper part shows Log(F (l)) computed using the DFA on the δ n , these plots have been displaced to avoid overlap. The lower part shows the Log(F (l)) from the DFA on the detrended δ ′ n . Plots for different degrees of the polynomial fit used to unfold the quantum spectra are shown in both cases. The dashed line shows the theoretical result for the GOE.

FIG. 6 :

 6 FIG. 6: Color online. Value of the α exponent of the power law behavior in the ensemble average calculation of the DFA for the fluctuation as a function of the degree of the polynomial fit used to unfold the quantum spectra. Two cases are shown, the α exponent values for the δ n fluctuation (filled squares, red online) and the case for the δ ′ n fluctuation obtained detrending the previous one (empty circles, blue online). The dashed line shows the theoretical result for the GOE. The error bars show the standard deviation of the ensemble.

VI. ACKNOWLEDGMENTS

This work was supported in part by grants from CONACyT, the Czech Science Foundation (202/09/0084) and the Czech Ministry of Education (MSM 0021620859).


José Mennesson

Christophe Saint-Jean

Laurent Mascarilla

NEW GEOMETRIC FOURIER DESCRIPTORS FOR COLOR IMAGE RECOGNITION

Keywords:

. In this paper, new generalized color Fourier descriptors (GCF D) are proposed. They depend on the choice of a bivector B acting as an analysis plane in a colorimetric space. The relevance of proposed descriptors is discussed on several color image databases. In particular, the influence of parameter B is studied regarding the type of images. It appears that proposed descriptors are more compact with a lower complexity and better classification rate.

INTRODUCTION

In the literature, there's many recent advances in terms of image recognition. The recognition process depends highly on discriminative and invariant descriptors. Among them, one can cite moment-based descriptors [START_REF] Flusser | Moments and Moment Invariants in Pattern Recognition[END_REF] such as Hu invariants, Legendre moments or Zernike moments. Beside this approach, SIFT (Scale-Invariant Feature Transform) descriptors are a popular choice giving very good results [START_REF] Lowe | Object recognition from local scale-invariant features[END_REF].

An alternative to these methods is to define descriptors in the frequency domain. In this framework, our paper concerns an extension of Fourier descriptors. Clearly, Fourier coefficients don't respect the classical invariances (translation, rotation and scale) and must be processed to obtain invariant descriptors. This paper propose an extension of a recent advance concerning generalized Fourier descriptors defined by F. Smach et al. [START_REF] Smach | Generalized fourier descriptors with applications to objects recognition in svm context[END_REF].

The extension of these descriptors to the color images is generally based on a marginal processing of the three channels R, G and B. Then, descriptors extracted from each channel are concatenated to form the description vector. In order to avoid this marginal processing, our proposal is to extract descriptors from a color Clifford-Fourier transform as defined by Batard et al. [START_REF] Batard | Clifford fourier transform for color image processing[END_REF]. This paper is structured as follows : In the section 2, the definition of the Clifford-Fourier transform is introduced. Then, section 3 presents the extraction of generalized descriptors from 2D Fourier transform. New geometric Fourier descriptors, corresponding to our contribution, are defined in section 4. Finally, experimentations on well-known databases are presented in the section 5.

A COLOR CLIFFORD-FOURIER TRANSFORM

The usual Fourier transform can be applied to L 2 (R m ;R or C) functions. Defining a color image as a function from R 2 to R 3 , the usual definition is no more usable. Several authors have proposed to tackle this problem by an ad-hoc processing of each channel. As an example, Smach et al. compute three Fourier transforms on each channel while Ell and Sangwine [START_REF] Ell | Hypercomplex fourier transforms of color images[END_REF] considered only two Fourier transforms: one on the luminance and the other on the chrominance plane. Intuitively, these works emphasize the notion of analysis direction in applications to color images.

Recently, Batard et al. [START_REF] Batard | Clifford fourier transform for color image processing[END_REF] have defined a Fourier transform for L 2 (R m ;R n ) functions. This one is mathematically rigorous and clarifies relations between the Fourier transform and the action of the translation group through a spinor group. They show that the previously proposed generalizations for color images (i.e. n=3) are particular cases of their definition.

In this paper, only the particular case m 2 and n 3 is considered and described shortly in the following.

Firstly, in the equation of the classical 2D Fourier transform, the term e ¡ipux vyq rotates f px, yq in the complex plane C. From a mathematical point of view, it corresponds to the action of the group S 1 on C which can be identified as the group action of SOp2q on R 2 . In order to generalize this principle to color images, one have to consider the action of the matrix group of SOp3q on R 3 . However, as described in [START_REF] Batard | Clifford fourier transform for color image processing[END_REF], a general and elegant expression may be written if the function corresponding to the image is embedded in the Clifford algebra R 4,0 .

Clifford algebra R 4,0 is a geometric algebra over the vector space R 4 , equipped with an euclidean quadratic form, on which a new product, called geometric product [START_REF] Hestenes | New algebraic tools for classical geometry[END_REF] is defined.

Then, a color image can be viewed as a R 2 into R 4,0 function: f px, yq rpx, yqe 1 gpx, yqe 2 bpx, yqe 3 0e 4 Within this framework, the rotation of a vector v by the angle θ, in the plane generated by a unitary bivector B can be written as

v Ñ e θ 2 B v e ¡ θ 2 B
Bivector B is an element of the geometric algebra R 4,0 which parametrizes the analysis direction. Function f can be decomposed as the sum of its parallel projection f B on the plane defined from B and its perpendicular projection f uB on the plane defined from I 4 B (I 4 being the pseudo-scalar of R 4,0 ). Skipping some technical details, one can prove that the values of f px, yq can be rotated by two independent rotations of the same angle, leading to the definition of a Clifford-Fourier transform:

x f B pu, vq » R 2 e pux vyq 2 B f B px, yq e ¡ pux vyq 2 B dxdy » R 2 e pux vyq 2 I4B f uB px, yq e ¡ pux vyq 2 I4B dxdy (1)
Identifying the plane defined from B and I 4 B to the complex plane C, this color Fourier transform can be computed easily and efficiently with two complex 2D-FFT:

x f B y f B y f uB

GENERALIZED FOURIER DESCRIPTORS

Generalized Fourier descriptors introduced by [START_REF] Smach | Generalized fourier descriptors with applications to objects recognition in svm context[END_REF] are defined from the group action of M 2 . This latter group is composed of translations and rotations on the plane. Two kinds of descriptors have been defined :

• "Spectral densities"-type invariants:

I r 1 pfq » 2π 0 § § § p f pr, θq § § § 2 dθ
• "Shift of phases"-type invariants:

I ξ1,ξ2 pfq » 2π 0 p f pR θ pξ 1 ξ 2 qq p f pR θ pξ 1 qq p f pR θ pξ 2 qqdθ
where f is the image, p f pr, θq is the Fourier transform expressed in polar coordinates in the frequency plane, ξ 1 and ξ 2 are variables of the frequency plane and R θ is a rotation of angle θ. It must be emphasized that, by construction, I r 1 and I ξ1,ξ2 are strictly invariant in R 2 with respect to the action of M 2 . Then, the descriptor vector for the first family of invariants, namely I r 1 , is defined as follows :

GF D1pf q 4 I 0 1 pfq § § § f p0, 0q § § § 2 , I 1
1 pfq I 0 1 pfq , . . . , I m 1 pfq I 0 1 pfq B
where m is the number of computed descriptors. In the same way, we define the descriptor vector corresponding to the second family of invariants I ξ1,ξ2 . These ones are called GF D2.

CONTRIBUTION

In order to deal with color images, a commonly used approach consists in computing descriptors on each color channel separately. Then, they are concatenated into a unique vector (e.g. [START_REF] Smach | Generalized fourier descriptors with applications to objects recognition in svm context[END_REF]). This method implies three FFT and three sets of descriptors. However, this marginal processing induces a loss of colorimetric informations that can be avoided by using color Clifford-Fourier transform.

Using the color Clifford-Fourier transform and the extraction of generalized Fourier descriptors, new color Fourier descriptors are defined as in Figure 1. The same construction based on I ξ1,ξ2 leads to GCF D2 B .

EXPERIMENTS

In order to evaluate the discrimination capacity of our proposed descriptors, the Coil-100 [START_REF] Nene | Columbia object image library (coil-100)[END_REF], the color FERET [START_REF] Phillips | The feret database and evaluation procedure for face recognition algorithms[END_REF] and the SFU [START_REF] Barnard | A data set for colour research[END_REF] databases are considered.

The COIL-100 database is composed of 7200 color pictures of size 128 ¢ 128 of 100 objects. Each picture has been taken with a black background and 72 different angles of view. This database, used in similar works [START_REF] Smach | Generalized fourier descriptors with applications to objects recognition in svm context[END_REF], can be qualified as "easy" from an image classification point of view.

The color FERET database is composed of face images of 1408 different persons, taken from different angles of view. In our tests, a set of 2992 images containing 272 persons equally represented by 11 pictures is selected and size of images is reduced to 128 ¢ 128. This database is more difficult than the first one.

The noisy color FERET database contains the precedent color FERET database, but a gaussian noise is added to the images. The parameter σ is fixed to 0.23, which is the maximum noise used in [START_REF] Smach | Generalized fourier descriptors with applications to objects recognition in svm context[END_REF].

The SFU database contains 20 objects photographed under 11 different illuminations and views. The size of images is also reduced to 128 ¢ 128. It is used to test the robustness of our descriptors under illumination variations.

Methodology

Size of descriptor vector. For the three databases, the descriptor vector length is set to 64 for each computed FFT. More precisely, GF D1 are built from 64 values of radius r in I r 1 and GF D2 are built from equally spaced values of ξ 2 in its polar domain r0, 2πs ¢ r1, 8s and ξ 1 set to p0, 1q. Then, after concatenation, this results into 192 descriptors in the marginal approach, and in 128 descriptors for the parallel plus orthogonal part.

Classification. The classification step is performed by the SVM (Support Vector Machine) algorithm as implemented in the LIBSVM library [START_REF] Chih | LIBSVM: a library for support vector machines[END_REF]. The parameters maximizing the recognition rates are σ 0.1 (gaussian kernel) and C 97.

To validate this decision step, a 10-fold cross-validation algorithm [START_REF] Jain | Data clustering: a review[END_REF] is performed. Moreover, a modified SFFS algorithm (Sequential Floating Forward Selection) [START_REF] Jain | Data clustering: a review[END_REF] is used to select bivectors B from their associated descriptor vectors.

Results

In our experiments, descriptors are built from unit bivectors B c e 4 , where c is a color to be chosen and the exterior The colors used to build the hundred random bivectors appears on each plot. Surprisingly, the rate obtained from B r B g B b is lower than that obtained from the best marginal B r . One can suppose that there are redundant data in the three descriptor vectors which reduce the classification capacity of the SVM. The great variability of the results proves that the recognition quality depends clearly on the choice of the bivector. Our experiments show that the bivector is hard to choose for a set of heterogeneous images. However, for the color FERET database, it appears that a good choice corresponds to the background color of the images. Indeed, it yields a high energy and discriminative Fourier coefficient for most of the images.

The different recognition rates computed from the descriptors previously defined on different databases are presented in the Table 1.

Firstly, the top of the table concerning results obtained from one bivector is studied. Immediately, one can see that the recognition rates permits to define an order on the descriptors such as GF D1 GCF D1 GF D2 GCF D2. Moreover, it can be noticed that on the noisy color FERET database, GF D2 descriptors give better results and are more robust against noise than GCF D1 descriptors with a half description length. Concerning the choice of the bivector, one can remark that B µ (luminance plus chrominance) is not a good choice for GCF D1. In the case of experimentations done from hundred random bivectors, the standard deviation obtained confirms that there is a relation between the recognitions rates and the choice of the bivector. The bottom of the table presents results achieved from a choice of a triplet of bivectors. The order previously defined is no more valid. However, the selection of a triplet of bivectors with an SFFS algorithm does not improve significantly results for GF D1 and GF D2 while results obtained from GCF D1 and GCF D2 are really better.

Using the SIFT implementation given by [START_REF] Vedaldi | VLFeat: An open and portable library of computer vision algorithms[END_REF][START_REF] Muja | Fast approximate nearest neighbors with automatic algorithm configuration[END_REF],we can compare our descriptors with SIFT descriptors. One can see that the results are lower than results achieved from GCFD1 or GCFD2 (Table 1). For the SFU database, it is due to the lack of interest points that can be extracted from these images.

CONCLUSION

We can conclude that we defined new color Fourier descriptors from a Clifford-Fourier transform [START_REF] Batard | Clifford fourier transform for color image processing[END_REF] and from generalized descriptors [START_REF] Smach | Generalized fourier descriptors with applications to objects recognition in svm context[END_REF]. We obtain much better recognition rates with GCF D1 and GCF D2 than those with the marginal approach for many choices of bivector B. Their computational costs are still low and the descriptor vector size is only twice larger. Some robustness against gaussian noise is obtained with GCF D2 descriptors. On chosen image databases, proposed descriptors gave higher recognition rates than the original SIFT descriptors. Unfortunately, selecting a good bivector remains still an open problem. However, x f B does not depend on bivector B. A deep work should be to defined GF D of Smach et al. directly for multi-vector valued functions. Finally, these first promising results encourage us to go farther with the general definition of the Clifford-Fourier transform [START_REF] Batard | Clifford fourier transform for color image processing[END_REF]. Indeed, this one is able to process 4D signals like for example color+infrared images.

( 2 )

 2 Note that a unit bivector B can be obtained from the geometric product of two unit orthogonal vectors as C 1 e 4 or C 2 C 3 (where C 1 ,C 2 and C 3 are RGB colors). These two choices of bivectors are equivalent because the orthogonal bivector to the unitary bivector B C 1 e 4 can be written as C 2 C 3 and vice versa. Obviously, an inverse transform exists and it is denoted by | f B [1].

Fig. 1 .f 2 dθ

 12 Fig. 1. Descriptors extraction from the color Clifford-Fourier transform.

Finally, the descriptor vector length is 2

 2 ¢ m: GCF D1 B pfq 2 GCF D1 B pfq, GCF D1 uB pfq @

product 1 Fig. 2 .

 12 Fig. 2. color FERET: Recognition rate with GCF D1 for 100 random bivectors

Table 1 .

 1 Recognition rates in % with GF D1, GCF D1, GF D2, GCF D2 and SIFT descriptors (last line)

		COIL-100	color FERET		Noisy color FERET		SFU
	Bivectors	GF D1	GCF D1	GF D1	GCF D1	GF D1	GCF D1	GF D2	GCF D2	GF D1	GCF D1
	Br	98.04	99.83	76.70	87.90	45.32	71.05	73.46	83.49	92.73	97.27
	Bv	98.06	99.56	73.66	79.65	46.83	61.99	75.26	78.64	93.18	93.64
	B b	96.90	99.86	70.49	84.49	48.49	73.46	74.77	81.78	95.00	95.45
	Bµ B rand p¢100q	98.49 98.42¨0.3	99.25 99.54¨0.3	73.03 73.72¨1.0	78.10 85.34¨2.9	55.28 54.23¨1.7	62.03 69.64¨3.2	77.34 76.59¨0.7	80.98 82.56¨1.8	95.45 94.56¨1.5	93.18 93.00¨1.6
	max.	98.87	99.89	76.14	90.37	57.55	77.27	78.41	87.00	97.27	96.36
	Br Bv B b	99.9	99.92	88.03	85.53	73.16	72.16	83.25	81.12	96.36	93.63
	B1 B2 B3(SFFS)	99.86	99.96	85.46	93.15	71.52	80.62	83.36	88.24	97.27	95.45
	SIFT	100		87.55		NA			91.82

When

vectors are orthogonal their geometric product reduces to their exterior product.


Guillaume Perillat-Merceroz

Patrice Gergaud

Pascal Marotel

Stephane Brochen

Pierre-Henri Jouneau

Guy Feuillet

Formation and annealing of dislocation loops induced by nitrogen implantation of ZnO

Keywords: numbers: 61.72.U-, 61.72.Dd, 68.55.Ln ZnO, ion implantation, annealing, semiconductor doping, transmission electron microscopy

Although zinc oxide is a promising material for the fabrication of short wavelength optoelectronic devices, p-type doping is a step that remains challenging for the realization of diodes. Out of equilibrium methods such as ion implantation are expected to dope ZnO successfully provided that the non-radiative defects introduced by implantation can be annealed out. In this study, ZnO substrates are implanted with nitrogen ions, and the extended defects induced by implantation are studied by transmission electron microscopy and X-ray diffraction (XRD), before and after annealing at 900 • C. Before annealing, these defects are identified to be dislocation loops lying either in basal planes in high N concentration regions, or in prismatic planes in low N concentration regions, together with linear dislocations. An uniaxial deformation of 0.4% along the c axis, caused by the predominant basal loops, is measured by XRD in the implanted layer. After annealing, prismatic loops disappear while the density of basal loops decreases and their diameter increases. Moreover, dislocation loops disappear completely from the sub-surface region. XRD measurements show a residual deformation of only 0.05% in the implanted and annealed layer. The fact that basal loops are favoured against prismatic ones at high N concentration or high temperature is attributed to a lower stacking fault energy in these conditions. The coalescence of loops and their disappearance in the sub-surface region are ascribed to point defect diffusion. Finally, the electrical and optical properties of nitrogen-implanted ZnO are correlated with the observed structural features.

I. INTRODUCTION

Because of its direct band gap of 3.37 eV, zinc oxide has gained a renewed interest during the last decade due to a wide range of possible applications in short wavelength optoelectronics. ZnO is an interesting alternative to gallium nitride (which has approximatively the same band gap) for several reasons. The high binding energy of the exciton (59 meV for ZnO versus 30 meV for GaN) could provide better luminescence at room temperature. Contrary to GaN-based devices which contain relatively rare materials like Ga and In, ZnO-based devices contain Zn and Mg which are abundant and non-toxic elements. Finally, relatively cheap ZnO substrates are commercially available contrary to GaN ones. This allows homoepitaxial ZnO layers to be grown with less structural defects than GaN heteroepitaxial layers. Important progresses have been made in the growth of ZnO and its alloys under the form of bulk crystals, thin layers, and nanowires. 1 But one of the remaining issues for the fabrication of ZnO-based devices is obtaining stable and reproducible p-type doping. 2 Group V elements such as N, P, As, or Sb substituting for O atoms are envisaged for p-type doping. Among them N is theoretically the best candidate a) Electronic mail: guillaume.perillat-merceroz@cea.fr because of its higher electronegativity and its ionic radius similar to that of O. However, from recent ab-initio calculations, the ionization energy could be 0.33 eV, 3 or 1.3 eV, 4 indicating that p-type doping with N may be difficult at thermodynamic equilibrium. Nevertheless, ZnO-based LEDs grown by molecular beam epitaxy at a temperature around 900 • C (a relatively low temperature compared to the melting point of 1975 • C) have been demonstrated using nitrogen for p-type doping. 5,6 Ion implantation is also promising for p-type doping of ZnO with nitrogen: it is a far from equilibrium method as a high supersaturation of point defects is created at room temperature, and N concentration can be beyond its equilibrium solubility. Moreover this method offers a precise doping through the control of dose, area localization, and dopant depth. For these reasons N + implantation of ZnO (called hereafter N implantation) is actively studied. [7][8][9][10][11] However implantation of energetic ions introduces defects in the crystal lattice which degrade the electrical 10 and optical 12 properties. Annealing is needed to remove these defects and to activate dopant atoms by moving them to the right lattice sites, keeping in mind that annealing at too high temperatures may on the contrary deactivate dopants by approaching equilibrium or by forming molecular nitrogen. 13 Although N-implanted and annealed ZnO remains n-type in general, 9,14,15 donoracceptor pair emission has been observed by photolu-minescence (PL). 8,9,16,[START_REF] Marotel | The 6th international workshop on zinc oxide an related materials[END_REF] Electroluminescence has been sometimes observed, but it is not clear if this was due to the formation of an actual p-n junction, 7,11 or to a metalinsulator-semiconductor junction. 10 These observations call for a thorough study of implantation defects and of their evolution upon annealing, in order to correlate the optical and electrical properties with the structural ones. The defects induced in ZnO by irradiation vary according to the type of irradiation (ion or electron), and according to the implanted element. Voids are produced by As implantation, [START_REF] Coleman | [END_REF] and clusters by Co or Mn implantation. 19 Although ion implantation in ZnO does not usually introduce any amorphous layer because of the large ionicity of its chemical bounds, 20 Si implantation can cause ZnO amorphization. 21 Ion bombardment thinning of heavily doped ZnO with low voltage Ar + results in the formation of interstitial dislocation loops lying in basal planes with a Burgers vector b = c/2 + p. 22 Electron irradiation of a ZnO thin foil creates interstitial dislocation loops lying in basal planes with b = c/2, and lying in prismatic planes with b = a. 23 It seems that dislocation loops are also formed during Au implantation, but the authors note that a more detailed transmission electron microscopy (TEM) study would be needed. 21 Defects induced by N implantation of ZnO, and sometimes their evolution with annealing, have also been studied. A Rutherford backscattering spectrometry study indicates that the absolute concentration of defects is much lower in ZnO than in GaN. 24 From X-ray diffraction (XRD), it is found that the crystalline quality is recovered after an annealing at 800 • C. 7 Positron annihilation spectroscopy studies reveal the formation of vacancy clusters, and explain their behaviour upon annealing. 14,15,25,26 Nevertheless, the crystallographic nature of the extended defects induced by N implantation of ZnO has never been determined, and the annealing effect has not been observed directly by TEM. In this study, N implantations were carried out in order to study the formation of structural defects and their recovery upon annealing. TEM and XRD were used to study an as-implanted sample, and a sample after implantation and annealing. Extended defects are identified to be dislocation loops lying in basal or prismatic planes, together with linear dislocations. The evolution of the defects with annealing was observed. Induced deformations in the implanted layer were measured by XRD and correlated to the observed microstructure. Considerations on dislocation energy and on local N concentration give some elements to discuss the different habit planes of the loops and the annealing effects. Optical and electrical characterizations obtained on the same samples are correlated with these structural observations. N being a potential p-type dopant for ZnO, the knowledge of the type of defects and their evolution with annealing is indeed of great importance in order to optimize structural, electrical, and optical properties of Nimplanted ZnO.

II. EXPERIMENTAL DETAILS

Two hydrothermally grown single-crystals of ZnO provided by Crystec were studied: one is an as-implanted substrate, and the other one is an implanted and annealed substrate. N + implantations were carried out at 0 • C on the oxygen face (i.e. the (0001) crystal face) with a NV-8200P Axcelis implanter. Samples were tilted 7 • relative to the ion beam direction to minimize channeling. The beam flux was 1.15 × 10 12 cm -2 s -1 . Three successive implantations at different energies and doses were carried out on each of the two samples in order to obtain a flat nitrogen profile. Accelerating voltages of 50, 120, and 200 kV with respective doses of 4 × 10 14 , 8 × 10 14 , and 10 15 cm -2 were used. The second sample was annealed at 900 • C during 15 minutes under O 2 at atmospheric pressure. These annealing conditions were chosen because they lead to the maximum intensity for the donor-acceptor pair emission in PL spectra. [START_REF] Marotel | The 6th international workshop on zinc oxide an related materials[END_REF] TEM sample preparation was carried out without ion milling in order to avoid artefacts, since Ar + milling of ZnO at 5 kV is known to introduce defects. 22 Although they can be minimized by a final milling at low voltage, 27 the wedge polishing technique was preferred. 28 TEM images were taken on a Jeol 4000 EX microscope operated at 400 kV. Scanning TEM images were taken on a FEI Titan microscope fitted with a probe C S -corrector and operated at 300 kV. XRD experiments were carried out on a Panalytical Xpert MRD diffractometer implementing triple axis measurements with the use of a fourfold Ge(220) monochromator, a threefold Ge(220) analyser, and Cu Kα radiation (1.540598 Å).

The following crystallographic notations are used in this paper for hexagonal ZnO: c = 0001 , a = 1 3 1120 , p = 1 3 0110 . {0001} planes are called basal planes. {1120} and {1100} planes are called prismatic planes. Dislocation loops lying in basal or prismatic planes are called basal or prismatic loops respectively. In the wurtzite structure, two types of basal loops can be found: b = c/2 dislocation loops which enclose extrinsic stacking faults with a stacking sequence ABABACBABA, and b = c/2 + p dislocation loops which enclose type one intrinsic stacking faults with a stacking sequence ABABABACACA. 22,[START_REF] Hirth | Theory of Dislocations[END_REF] III. RESULTS

A. N concentration profiles

Fig. 1 shows the implantation profiles measured by secondary ion mass spectrometry (SIMS) before and after annealing, together with the simulated profile obtained with the SRIM software. [START_REF] Ziegler | [END_REF] SIMS measurements and the simulation are in good agreement.The annealed and the as-implanted sample present nearly the same concentration profile. Thus N diffusion from the implanted zone in FIG. 1. N concentration profiles simulated with the SRIM software, [START_REF] Ziegler | [END_REF] and measured by SIMS on the as-implanted sample and on the implanted and annealed sample.

ZnO or to the surface appears to be negligible at 900 • C in our experimental conditions.

B. Crystallographic defects induced by N implantation

Fig. 2 shows cross-sectional TEM images of defects in the as-implanted layer. Three types of defects can be identified: dislocation loops lying in basal planes; dislocation loops lying in prismatic planes, and linear dislocations.

Basal loops

On the [0002] weak-beam image [Fig. 2(a)], overlapping black and white contrasts due to a high density of basal loops are visible at depths between 0 and 400 nm. These depths correspond to the zone of high N concentration [Fig. 1]. The habit plane of the loops in this zone was identified on cross sectional high resolution images taken along a 1120 zone axis, as shown in Fig. 2(b). An additional (0001) insertion plane is detectable, indicating that the loops lie in basal planes and are of interstitial type. The mean loop diameter was determined to be 5 nm from various high resolution images. Due to the high density and the small size of these basal loops, their Burgers vector is difficult to determine without ambiguity from a g.b 31 or in materials irradiated by electrons, neutrons or ions such as Zr, 32 Co, 33 SiC, 34 and ZnO. D. Deformations in the implanted zone as measured by XRD Fig. 7 shows the results obtained by high resolution XRD. Fig. 7 (a) displays θ/2θ scans done on the symmetrical [0002] reflection for an unimplanted sample, the asimplanted sample, and the implanted and annealed sample. Fringes on the left of the [0002] peak indicates that nitrogen implantation in ZnO induces a deformed layer with a larger c parameter. Reciprocal space maps on the asymmetrical [1124] reflection [Fig. 7 (b)] demonstrate that the dilatation of the implanted layer occurs only in the c direction, with no in-plane deformation. This gives evidence of an increase of the ZnO unit cell induced by implantation. Goland and Keating have already demonstrated that interstitial dislocation loops in hexagonal structures yield an increase of the lattice parameter in the direction perpendicular to the habit plane of the loops. 35 This phenomenon, intensively studied in nuclear material science, was called "irradiation growth". 36 In our study, we infer that the basal loops yield an increase of the c parameter. The deformation induced by the prismatic loops is not detectable, as deduced from Fig. 7 (b). Simulations of rocking curves based on kinematical diffraction are shown in Fig. 8(c). The experimental rocking curve is initially fitted with a uniform deformation profile as represented in the inset of Fig. 8(a). This corresponds to a layer of 400 nm with a 0.4% uniaxial deformation. How- ever, a better fit of the experimental rocking curve is obtained with a deformation profile which follows qualitatively the implantation profile, as represented in the inset of Fig. 8(b). This suggests that the distribution of the basal loops follows roughly the implantation profile. The XRD scan for the annealed sample is also presented in Fig. 7(a): a small deformation of only 0.05% remains after annealing, which can be ascribed to the remaining basal loops.

IV. DISCUSSION

Table I sums up the TEM observations of implantation defects obtained by Yoshiie et al. 23,37 for electron irradiation of pure ZnO thin foils, by Courderc et al. 22 for Ar implantation of Bi-Mn-Ti-doped ZnO, and by ourselves for N implantation. All the observed basal dislocation loops are of interstitial type, whereas prismatic loops are only assumed to be interstitial.

A. Formation of dislocations during implantation 1. Formation of interstitial dislocation loops

First-principles calculations show that native point defects are very mobile in ZnO. [38][39][40] Calculated migration barriers and the corresponding temperatures at which native defects become mobile (considering that a point defect is mobile if its jump rate is more than 1 per second) are reported in Table II. 38 It appears that at room temperature, interstitials are very mobile, and vacancies almost immobile. In ion implanted ZnO, Frenkel pairs are formed by a cascade damage. Some of them annihilate, but mobile interstitials can agglomerate to form dislocation loops. Immobile vacancies stay isolated or may form vacancy clusters of a few vacancies. 14,15,25,26 Our observations of interstitial loops are thus coherent with theoretical predictions. Concerning previous TEM investigations (reported in Table I), Yoshiie et al. and Couderc et al. also observed that basal loops were of interstitial type.

Type of dislocation loops according to N concentration

Let us discuss about the parameters that may determine the habit planes of interstitial loops. Föll et al. proposed a rule to predict the habit planes of dislocation loops in hexagonal closed-packed structures depending on the c/a ratio. 33 If d {0001} /d {1010} < 1 (where d is the lattice plane spacing), which is equivalent to c/a < √ 3 = 1.732, the loops should form in the {1010} prismatic planes. If c/a > 1.732, the loops should form in basal planes. Griffiths discussed more precisely the habit planes of interstitial loops in a review on the microstructure of irradiated hexagonal compact metals. 41 Föll's rule is respected for pure materials. But for impure materials with c/a < 1.732, both prismatic and basal loops are observed. In particular, dislocation loops formed by quenching in Mg (for which c/a = 1.623) were studied. 31 Prismatic loops where observed for Mg with 0.0001% impurities, while basal loops were observed for Mg with 0.05% impurities. In the ZnO case, prismatic loops should form preferentially, as c/a = 1.602. However, taking into account Griffiths' observations, the number of basal loops should increase relatively to the number of prismatic loops when the impurity concentration increases. This is in agreement with the results for irradiated ZnO summarized in Table I: prismatic loops together with basal loops are present in relatively pure material (Yoshiie et al.) whereas only basal loops are present in impure materials (Couderc et al.). In our study, basal loops are observed in the high N concentration region and prismatic loops in the low N concentration region, in agreement with Griffiths' observations. This dependence of the habit plane with the impurity concentration can be explained by considering the energy of a dislocation loop. The total elastic energy of a loop is E = E l × l + γ × S where E l is the energy per unit-length of the dislocation line, l the length of the dislocation, γ the energy of the stacking fault, and S its surface. A general result for a large variety of materials is that the stacking fault energy is lowered when the impurity concentration is increased. 42 For example, in the case of wurtzite materials like GaN, the density of basal stacking faults is increased by a factor of 10 with Si doping. 43 Ab-initio calculations indicate this is due to a lower basal stacking fault energy resulting from the Si doping. 44 As observed for many crystals (for example Si), 45 and for ZnO in particular, 23 loop diameters increase upon annealing while their density decreases. As two dislocation loops attract each other, loops tend to coalesce. [START_REF] Hirth | Theory of Dislocations[END_REF] Because the Burgers vector of the loops do not lie in the loop plane, this coalescence does not occur by glide. It occurs rather by climb, which is a diffusion-assisted phenomenon. The disappearance of all loops in the subsurface region (up to 65 nm under the surface) with annealing is observed. One could attribute this to the weaker density of dislocation loops in this region, as deduced from the XRD study (see III D). But it is not sufficient to explain the complete disappearance of the loops. The specific role of surface must be considered. Actually, it was already observed that the density of loops produced by implantation in silicon decreases when the implantation depth approaches the surface. 46 This was ascribed to the diffusion of interstitials towards the surface. Sub-surface area free of dislocation loops were also observed in hexagonal materials such as Zr alloy, 47 Mg, 31 and CdTe. 48 Because hexagonal materials are anisotropic, it was pointed out that all surfaces do not play the same role. 47 In ZnO, irradiated thin foils 23 and Ga-implanted ZnO nanorods 49 become free of dislocations after annealing at 800 • C. These examples stress the role of surfaces to remove dislocation loops. Surfaces act as a sink for point defects and thus causes the "dissolution" of dislocation loops. In this study, the disappearance of prismatic loops is observed for an annealing at 900 • C. Yoshiie et al. observed that the basal loops to prismatic loops ratio increased with the irradiation temperature from 200 • C to 500 • C. Extrapolation of their results shows that no prismatic loops should be present after 720 • C, which is coherent with our study. Prefer-ential nucleation of basal loops was also observed in Zr alloys at higher temperature. 32,41 A lower basal stacking fault energy due to a higher temperature 42 may cause prismatic loops to be unstable and disappear.

Considerations to improve structural recovery

After annealing, the XRD analysis indicates a good structural recovery with a residual deformation of only 0.05%. However, PL spectra show no donor-acceptor pair emission when carried out directly on the sample surface. But when the top-surface layer (those without defects on a depth of 65 nm) is removed by chemical-mechanical polishing, donor-acceptor emission is evidenced by PL. [START_REF] Marotel | The 6th international workshop on zinc oxide an related materials[END_REF] Hall measurement using the Van der Paw configuration was realized on the annealed sample after removal of the surface layer. However, the implanted layer was found too resistive to be measured, most probably because of the remaining dislocations observed in TEM. In our view, the electroluminescence observed in the literature for Nimplanted ZnO annealed at 800 • C is possibly due to a metal-insulator-semiconductor junction. 7 To conclude, a better structural recovery is still needed to produce lightemitting devices with N implanted ZnO as a p-type layer. Several solutions may be envisaged for a better structural recovery of the implantation induced defects.

• A longer annealing, or an annealing at higher temperature are possible, but it would lead to the disappearance of donor-acceptor pair emission, [START_REF] Marotel | The 6th international workshop on zinc oxide an related materials[END_REF] and to the possible deactivation of N doping by forming N 2 molecules. 13 • Lower implantation doses would produce less defects, and consequently the annealed material would have less defects too. Moreover, more prismatic loops should form for low N dose, and they should disappear with annealing. But the dose must remain sufficient to incorporate enough active dopants.

• Finally, implantation could be carried out on ZnO nanowires to take advantage of a better recovery due to the larger surface to volume ratio in these structures. Actually it was observed that structural defects were removed in Ga implanted nanowires after an annealing at 800 • C. 49 Moreover, UV electroluminescence was obtained from a p-n junction in ZnO nanowires, with arsenic implantation doping. 50 V. SUMMARY N implantation in ZnO at room temperature produces immobile vacancies and mobile interstitials. Interstitials agglomerate into dislocation loops lying in basal and prismatic planes. For depths between 0 and 400 nm, where N concentrations are high, loops lying in basal planes are observed. Deeper than 400 nm, where N concentration decreases, loops lying in prismatic planes with b = a are observed. Linear dislocations with b = a are also observed. The implanted layer is shown by XRD to have a larger c parameter, which is ascribed to the predominant basal loops. After annealing at 900 • C, basal loops with b = c/2 + p and linear dislocations with b = a are still observed, but prismatic loops disappear. A mechanism is proposed to explain these phenomena: basal loops are favoured against prismatic loops at high N concentration or high temperature because the basal loop stacking fault energy decreases. With annealing, basal loops coalesce: they are less numerous and bigger. They disappear below the surface. Both phenomena are ascribed to diffusion of point defects. Although XRD and PL results seems to show a good recovery of the material, the remaining basal loops observed in TEM lead to a non conductive layer, as demonstrated by Hall effect measurements. Propositions are made to improve the structural recovery of N-implanted ZnO.

FIG. 2 .

 2 FIG. 2. TEM weak-beam images of the as-implanted sample observed in cross-section with (a) g = [0002], and (c) g = [1100] diffracting conditions. (b) High resolution TEM image taken along the 1120 zone axis in the area indicated in (a).

  analysis or from high resolution images. On the [1100] image [Fig.2(c)], the basal loops show a faint contrast. It may be a residual contrast (g • (b × u) = 0 even if g.b = 0, where u is a unit vector along the dislocation line), or it may correspond to a situation where g.b = 0. This is compatible with a Burgers vector b equal either to c/2 or to c/2 + p, which are the two possibilities for interstitial loops lying in basal planes in the hexagonal structure.

FIG. 3 .

 3 FIG. 3. Cross sectional weak-beam images of the implanted zone at 400 nm under the surface with (a) g = [0002], and (b) g = [1100] diffracting conditions.

23

 23

3 .

 3 FIG. 5. Weak-beam images of N implanted ZnO after annealing with (a) g = [0002], (b) g = [1010], and (c) g = [2110] diffracting conditions.

[

  Fig. 2(a)]. Consequently their Burgers vectors lie in the basal plane. They are also visible on plane-view images taken near the 0001 zone axis [Fig. 4]. Two linear dislocations are visible for g = [0110] and g = [1100], and are out of contrast for g = [1010]. It results from this analysis that the Burgers vector of these linear dislocations is b = 1/3[1210] = a. The mechanism for the formation of these linear dislocations is not clearly understood yet.

FIG. 6 .

 6 FIG. 6. (a) Low magnification bright-field scanning TEM image of basal loops after annealing, and (b) high-resolution scanning TEM image of the side of one loop (corresponding to the area marked by a rectangle in (a)) taken along a 1120 zone axis with a high-angle annular dark-field detector.

FIG. 7 .

 7 FIG. 7. XRD results. (a) θ/2θ scans on the symmetrical [0002] reflection obtained on ZnO substrates before implantation, after N implantation, and after implantation and annealing. (b) Reciprocal space map on a [1124] asymmetric peak obtained on the as-implanted sample

FIG. 8 .

 8 FIG. 8. (c) θ/2θ scan on the symmetrical [0002] reflection obtained on the as-implanted ZnO substrate, and XRD simulations for various d-spacing gradients along the c axis: (a) uniform deformation from 0 to 400 nm, and (b) deformation profile proportional to the implantation profile (see Fig. 1).

b = c/ 2 b

 2 = c/2 + p b = c/2 or b = c/2 + p (ambiguity) for high CN region Prismatic loops with b = a Prismatic loops with b = a for low CN region Implanted and annealed sample Fewer prismatic loops. No more loops if T > 800 • C. Interstitial basal loops with b = c/2 + p. Prismatic loops and loops close to the surface are removed at 900 • C.

  Because prismatic loops have a Burgers vector b = a, and thus no stacking fault inside (the dislocation is perfect), their energy does not change when the impurity concentration increases. But the basal loops with b = c/2 + p or b = c/2 contain a stacking fault: their energy is lowered when the impurity concentration is increased. Thus basal loops become more stable than prismatic loops for high impurity concentration. The Burgers vectors of basal loops may also be accounted for by the impurity concentration. In our study, and in Couderc et al.'s one with impure ZnO, basal loops have b = c/2 + p with a type one intrinsic stacking fault. In Yoshiie et al.'s study with pure ZnO, the Burgers vector is b = c/2 with an extrinsic stacking fault. To explain these experimental observations, we propose that a high impurity concentration may lower the intrinsic stacking fault energy more than the extrinsic one and thus favours b = c/2 + p in impure materials.B. Recovery of implanted ZnO with annealing 1. Behaviour of dislocation loops

TABLE I .

 I Types of dislocation loops observed by TEM in ion or electron irradiated ZnO before and after annealing. The N concentration is noted CN .

		Yoshiie et al. 23,37	Couderc et al. 22	This study
		Electron irradiation	Ar implantation	N implantation
		700 keV	5 keV	50 keV to 200 keV
		Pure ZnO	Zn97.4Bi1Mn1Ti0.6O	ZnO with 0 < CN < 0.22%
	Implanted	Interstitial basal loops with	Interstitial basal loops with	Interstitial basal loops with
	sample

TABLE II .

 II Migration barriers and temperature at which ZnO native defects become mobile (from Janotti et al.38 ).

	Zinc interstitial	0.57 eV (-54 • C)
	Oxygen interstitial	0.87 eV (62 • C) or 1.14 eV (166 • C)
	Zinc vacancy	1.40 eV (266 • C)
	Oxygen vacancy	2.36 eV (636 • C)

ACKNOWLEDGMENTS

The authors would like to acknowledge Carole Granier for her precious technical assistance, Frédéric Milesi for implantations, Jean-Paul Barnes for making some corrections to the English of this publication, and the french national research agency (ANR) for funding through the Carnot program (2006/2009).


Thanh Phuong Nguyen
email: nguyentp@loria.fr

Bertrand Kerautret
email: kerautre@loria.fr

Ellipse detection through decomposition of circular arcs and line segments

. The first step of such a detection is to locate ellipse candidate with a grouping process exploiting geometric properties of adjacent arcs and lines. Then, for each ellipse candidate we extract a compact and significant representation defined from the segment and arc extremities together with the arc middle points. This representation allows then a fast ellipse detection by using a simple least square technique. Finally some first comparisons with other robust approaches are proposed.

Introduction

Shape identification is an important task in image analysis. Ellipse is a basic shape that can appear naturally in images from 3D environment. Therefore, ellipse detection is a key problem in many applications in the field of computer vision or pattern recognition.

In general, we can group existing methods into three main approaches. The first one relies on the Hough transform [START_REF] Huang | Elliptical feature extraction via an improved hough transform[END_REF][START_REF] Aguado | On using directional information for parameter space decomposition in ellipse detection[END_REF][START_REF] Daul | From the Hough Transform to a New Approach for the Detection and Approximation of Elliptical Arcs[END_REF][START_REF] Chia | Ellipse detection with hough transform in one dimensional parametric space[END_REF][START_REF] Lu | Detection of incomplete ellipse in images with strong noise by iterative randomized Hough transform (IRHT)[END_REF]. These methods transform image into parametric space and then take the peaks in this space as candidate of ellipse. Generally, it requires a parameter space that has five dimensions -contrariwise to two for straight line detection and three for circle detection, so it needs more execution time and memory space than the two last approaches. Some modifications [START_REF] Daul | From the Hough Transform to a New Approach for the Detection and Approximation of Elliptical Arcs[END_REF][START_REF] Chia | Ellipse detection with hough transform in one dimensional parametric space[END_REF][START_REF] Lu | Detection of incomplete ellipse in images with strong noise by iterative randomized Hough transform (IRHT)[END_REF] of Hough transform have been proposed to minimize storage space and computation complexity. Daul et al. [START_REF] Daul | From the Hough Transform to a New Approach for the Detection and Approximation of Elliptical Arcs[END_REF] reduce the problem to two dimensional parametric space. Later, Chia et al. [START_REF] Chia | Ellipse detection with hough transform in one dimensional parametric space[END_REF] introduced a method based on Hough transform in one dimensional parametric space. Lu et al. [START_REF] Lu | Detection of incomplete ellipse in images with strong noise by iterative randomized Hough transform (IRHT)[END_REF] proposed an iterative randomized Hough transform (IRHT) for ellipse detection with strong noise.

The second one uses least square fitting technique [START_REF] Ahn | Geometric least squares fitting of circle and ellipse[END_REF][START_REF] Ciobanu | Contour profiling by dynamic ellipse fitting[END_REF][START_REF] Fitzgibbon | Direct least square fitting of ellipses[END_REF][START_REF] Ahn | Least-squares orthogonal distances fitting of circle, sphere, ellipse, hyperbola, and parabola[END_REF] that minimizes the sum of square error. There are two main types of least square fitting (see [START_REF] Ahn | Least-squares orthogonal distances fitting of circle, sphere, ellipse, hyperbola, and parabola[END_REF]) that are characterized by the definition of error distances: algebraic fitting and geometric fitting. Concerning the first type, the error distance is defined by considering the deviation at each point to the expected ellipse described by implicit equation F (x, a) = 0 where a is vector of parameters. Contrariwise, for the second type, the error distance is defined as orthogonal distance from each point to the fitting ellipse.

The third group of approach detects ellipse candidates by using their moment [START_REF] Voss | Invariant fitting of planar objects by primitives[END_REF][START_REF] Lee | Moment preserving detection of elliptical shapes in gray-scale images[END_REF][START_REF] Heikkila | Moment and curvature preserving technique for accurate ellipse boundary detection[END_REF][START_REF] Zunic | Efficiency of characterizing ellipses and ellipsoids by discrete moments[END_REF][START_REF] Rosin | Measuring shape: ellipticity, rectangularity, and triangularity[END_REF].

We propose a new method for ellipse detection based on the decomposition of an edge image into arc and line primitives. The main contribution of this paper is to propose a pre-processing step that allows to speed up the detection of ellipse based on a linear scanning process on the sequence of arc and line primitives. The rest of this paper is organized as follows. The following section recalls a method for the representation of a digital curve by arcs and line segments. Section 3 presents the proposed method for ellipse detection before experimentation.

Descriptor based on arc and line primitives

In this section, we recall a linear method [START_REF] Nguyen | Decomposition of a curve into arcs and line segments based on dominant point detection[END_REF] for the decomposition of a digital curve into circular arcs and line segments.

Tangent space representation and properties of arc in the tangent space

Nguyen and Debled-Rennesson proposed in [START_REF] Nguyen | A linear method for segmentation of digital arcs[END_REF] some properties of arcs in tangent space representation that are inspired from Latecki [START_REF] Latecki | Shape similarity measure based on correspondence of visual parts[END_REF].

Let C = {C i } n i=0 be a polygon, l i -length of segment C i C i+1 and α i = ∠( ----→ C i-1 C i , ----→ C i C i+1 ). If C i+1 is on the right of ----→ C i-1 C i then α i > 0, otherwise α i < 0 (see illustration of Fig. 1(a)).
Let us consider the transformation that associates a polygon C of Z 2 to a polygon of R 2 constituted by segments T i2 T (i+1)1 , T (i+1)1 T (i+1)2 , 0 ≤ i < n (see Fig. 1(b)) with:

T 02 = (0, 0), T i1 = (T (i-1)2 .x + l i-1 , T (i-1)2 .y), i from 1 to n, T i2 = (T i1 .x, T i1 .y + α i ), i from 1 to n -1. Nguyen et al. also proposed in [START_REF] Nguyen | A linear method for segmentation of digital arcs[END_REF] some properties of a set of sequential chords of a circle in the tangent space. They are resumed by proposition 1 (see also Fig. 2).

α 1 α 2 α 3 C 0 C 1 C 2 C 3 C 4 (
1. [16] Let C = {C i } n i=0 be a polygon, α i = ∠( ----→ C i-1 C i , ----→ C i C i+1 ) such that α i ≤ α max ≤ π 4 . The length of C i C i+1 is l i , for i ∈ {1, .
. . , n}. We consider the polygon T (C), that corresponds to its representation in the modified tangent space, constituted by the segments

T i2 T (i+1)1 , T (i+1)1 T (i+1)2 for i from 0 to n -1. M pC = {M i } n-1 i=0 is the midpoint set of {T i2 T (i+1)1 } n-1 i=0 . So, C is a polygon whose vertices are on a real arc only if M pC = {M i } n-1
i=0 is a set of quasi collinear points.

From now on, M pC is called the midpoint curve.

O C 2 C 3 H 0 H 1 C 0 C 5 C 4 C 1 α 1 H (a)
A set of sequential chords of an arc.

x y 0

T (i+2)1 M i+1 T (i+1)2 T i2 M i T (i+1) 1
T (i-1)2 M i-1 T i1
(b) Its property in tangent space representation.

Fig. 2. The chords in tangent space.

Arc line decomposition

Proposition 1 can be used to decide if a digital curve is a circular circle by detecting straight line segment in the tangent space. Moreover, it is also used for the decomposition of a curve into arcs and line segments. Nguyen introduced the definition below.

Definition 1. In the curve of midpoints in the tangent space, an isolated point is a midpoint satisfying that the differences of ordinate values between it and one of its 2 neighboring midpoints on this curve is higher than the threshold α max . If this condition is satisfied with all 2 neighboring midpoints, it is called a full isolated point

Let us consider Fig. 3. In this example, there are all basic configurations among the primitive arc and line: arc-arc, arc-line and line-line. Fig. 4 presents these configurations in detail in the tangent space. Concerning the midpoint curve (M pC) in the tangent space, Nguyen et al. [START_REF] Nguyen | Decomposition of a curve into arcs and line segments based on dominant point detection[END_REF] introduced several remarks below.

-An isolated point in M pC corresponds to an extremity among two adjacent primitives in C. -A full isolated point in M pC corresponds to an line segment in C.

-An isolated point in M pC can be co-linear with a set of co-linear points that corresponds to an arc. In [START_REF] Nguyen | Decomposition of a curve into arcs and line segments based on dominant point detection[END_REF], Nguyen et al. proposed an algorithm (see3 algo. 1) to decompose a curve C into arcs and straight line segments. First, the sequence of dominant points (DpC) of C is computed by using an algorithm presented in [START_REF] Nguyen | A discrete geometry approach for dominant point detection[END_REF]. DpC is then transform in the tangent space and the M pC curve is constructed. An incremental process is then used and each point of M pC is tested: if it is not an isolated point (in this case, it corresponds to an arc segment in C), the blurred segment recognition algorithm [START_REF] Debled-Rennesson | Optimal blurred segments decomposition of noisy shapes in linear time[END_REF] permits to test if it can be added to the current blurred segment (which corresponds to an arc in C). If it is not possible, a new blurred segment starts with this point.

Ellipse detection

We present hereafter a new method for ellipse extraction from edge map of an image. It is based on three steps:

-Construction of a representation based on arc and line primitive of the edge map of input image. The first step is done by applying the decomposition of a curve into arcs and lines presented in the above section. We construct the corresponding edge image from input image by using Canny filter. This edge image is considered as a list of digital curves. Thanks to above technique [START_REF] Nguyen | Decomposition of a curve into arcs and line segments based on dominant point detection[END_REF], we can obtain a compact representation of this edge image based on arc and line primitives. Grouping of arcs and lines based on geometric property. In this step, ellipse candidate is detected by grouping adjacent arcs and lines based on geometric property. Its main idea is to group the adjacent primitives that have a same tangent vector at the common extremity. In practice it is done by verifying the angle between two adjacent primitives.

Let us consider two adjacent primitives. There are three possible configurations: arc-arc, arc-line and line-line. We define the angle among two primitives that depends on its configuration as follows.

-Arc-line: Let us see Fig. Thanks to the notion of angle among two adjacent primitives and sequential property of the representation based on arcs and line segments, a linear scanning process is used for grouping adjacent primitives satisfying that the angle among two adjacent primitives does not exceed a fixed threshold. A such group of primitives that contains at least one primitive of arc is called an ellipse candidate.

To avoid false positive with small detected ellipse, we use two constraints about the arcs in each group of primitives: the maximal radius of arc must be greater than 5 and the total subtending angle of arcs must be greater than π 5 .

Fitting of ellipse candidate. For each ellipse candidate constructed from the extracted arcs and lines, we try to fit it by using least square fitting. Contrariwise to existing techniques based on least square fitting, we need a very small set of extracted points for fitting. This good property is given by the repre-

A B C D E M 1 M 2
Fig. 6. Ellipse fitting based on least square fitting.

sentation of the curve by arc and line primitives. The set of extracted points for ellipse fitting is constituted from extremities of arcs and lines. Moreover, it contains also midpoints of arcs to reduce the approximated error between the fitting ellipse and the curve reconstructed by arcs and lines. Concerning least square fitting, there are two main categories: algebraic fitting and geometric fitting (see [START_REF] Ahn | Least-squares orthogonal distances fitting of circle, sphere, ellipse, hyperbola, and parabola[END_REF]). Let us consider a general conic described by this implicit function: F (A, X) = A•X = ax 2 +bxy +cy 2 +dx+ey +f = 0 where We have experimented the proposed method on a 2.26 GHz CPU linux computer, with a 4Go of RAM. The results are illustrated in Fig. 7 and Fig. 8. From the input image (a), the edge image (b) is computed. Then, the detected arcs and lines from algorithm 1 are represented in (c). Afterwards, a scanning process is applied to group arc and line primitives to detect the ellipse candidate (d) by verifying the angle between adjacent primitives. Finally, for each ellipse candidate, the technique of least square fitting is applied to construct the fitting conic on the small set of extracted points (marked points in detected ellipses in figures Fig. 7.e and Fig. 8.e). Fig. 7.f and Fig. 8.f present the results obtained by other methods [START_REF] Wu | Robust real-time ellipse detection by direct least-square-fitting[END_REF], [22].

Table 1 shows some information about processing of ellipse detection from the previous experiments. Thanks to the representation based on arcs and lines, the number of primitives for processing is reduced from 13699 (resp. 7778) to 308 (resp. 261) for image in Fig. 7 (resp. Fig. 8). In addition, the scanning process on the sequence of primitives is done in linear time due to its sequential property. After applying this process, the number of ellipse candidates is reduced to 27 (resp. 6) for image in Fig. 7 (resp. Fig. 8). Moreover, the average number of extracted points for ellipse fitting of each candidate is only 11.18 and 8.17.

The proposed method has two main advantages that guarantees its fastness:

-An efficient pre-processing step for the detection of ellipse candidates. It is based on a fast method to represent the edge image by arcs and line segments and a linear scanning process for the detection of ellipse candidates. -A small set of extracted points for ellipse fitting. Thanks to the representation based on arcs and lines, we don't need all of points corresponding to ellipse candidate for ellipse fitting. We only use extremity points of the primitives (arcs, lines) and midpoints of arcs for this task.

Conclusions

We have presented a promising new method for ellipse detection in images.

The scanning process for detection of ellipse candidate is efficient because the representation based on arcs and lines allows us to work with a small number of primitives in relation with the number of points in edge image. In addition, this process is done in linear time. Moreover, for each ellipse candidate, the least square fitting is not directly applied on all the points of the candidate. In future work we plane to perform other comparisons with different approaches and include a recent noise detection defined through the concept of meaningful scale [START_REF] Kerautret | Multi-scale analysis of discrete contours for unsupervised noise detection[END_REF].

Fig. 1 .

 1 Fig. 1. Tangent space representation

Fig. 3 .Fig. 4 .

 34 Fig. 3. An example of curve

Algorithm 1 : [ 1 ]-

 11 Curve decomposition into arcs and lines 3 . Data: C = {C1, . . . , Cn}-a digital curve, αmax-maximal angle, ν-width of blurred segments Result: ARCs-set of arcs, LIN Es-set of lines begin Use [18] to detect the set of dominant points: DpC = {D0, . . . , Dm}; BS = ∅; Transform DpC in the tangent space as T (DpC); Construct the midpoint curve M pC = {Mi} m-1 i=0 of horizontal segments of T (DpC); for i=0 to m-1 do {Ci} e i b i -part of C wich corresponds to Mi; if (BS ∪ Mi is a blurred segment of width ν [19]) and (|Mi.y -Mi-1.y| < αmax) then BS = BS ∪ Mi; else C ′ -part of C corresponding to BS; Push C ′ to ARCs; BS = {Mi}; if (|Mi.y -Mi+1.y| > αmax) then Push C b i Ce i to LIN Es; BS = {∅} ; end Grouping of arcs and lines for detection of ellipse candidate based on geometric properties. -Fitting of ellipse candidate based on least square fitting.

Fig. 5 .

 5 Fig. 5. Arc and line grouping based on geometric property.

  5.a. The arc has its center O 1 . Two primitives share a common point T. T z is perpendicular with T B. The angle among two primitives is define as ∠O 1 T z. -Arc-arc: Let us see figure 5.b. Two adjacent arcs whose centers are O 1 and O 2 share common point T. The angle among these arcs is define as ∠O 1 T O 2 . -Line-line: Let us see figure 5.c. The angle between two line segments AB and BC is π -∠ABC.

Fig. 6

 6 illustrates this strategy. The ellipse candidate is composed of segments AB, BC and arcs CD, DE. The data set contains A, B, C, D, E the extremities of arc and line primitives, M 1 , M 2 the midpoints of arcs CD and DE.

Fig. 7 .

 7 Fig. 7. Comparison with Wu's method [21].

Table 1 .

 1 Ellipse detection in images in figures 7, 8.

	Figure Image size N o of points	N o	of	N o of can-	N o of el-	N o of points	Times
			in edge image	primitives	didates	lipses	per candidate	(ms)
	7	584x440	13699	308		27	23	11.18	880
	8	508x417	7778	261		6	5	8.17	520

Note this algorithm includes some corrections of algo. 3 of[START_REF] Nguyen | Decomposition of a curve into arcs and line segments based on dominant point detection[END_REF].

By default, αmax = π 4 , ν = 1, βmax = π

(see algo. 2).

, xy, y 2 , x, y] t . F (A, X i ) is defined as algebraic (resp. geometric) distance to the conic F (A, X) = 0. The least square fitting is used to minimize the sum of squared error distance: m i=0 F (A, X i ) 2 with none trivial solution A = [0, 0, 0, 0, 0, 0] t . Many works have been proposed for minimizing this sum of square error. In our work, the ellipse is fitted from the extracted set of points by using the code of ellipse fitting proposed by Ohad Gal [START_REF] Gal | Matlab code for ellipse fitting[END_REF].

Proposed Algorithm. Algorithm 2 describes the ellipse detection based on the representation of an edge image by arc and line primitives. The scanning process is applied for the detection of ellipse candidates by grouping adjacent arcs and lines which have their adjacent angles smaller than a threshold (β max ). For each ellipse candidate, the algorithm tries to construct an approximated conic by using least square fitting on a small set of points constructed from extremities of arcs and lines and middle points of arcs. If the conic is an ellipse, it is considered as a detected ellipse. In practice, the threshold of angle is set from π 8 to π 5 .