Ferréol Soulez
email: ferreol.soulez@univ-lyon1.fr

Charles Crespy

Valerie Kaftandjian
email: valerie.kaftandjian@insa-lyon.fr

Philippe Duvauchelle

Valérie Kaftandjian

Diffraction peaks restoration and extraction in energy dispersive X-ray diffraction

Keywords: Energy Dispersive X ray Diffraction, Non Destructive Testing, spectrum restoration, inverses problems

come

Introduction

X-ray diffraction (XRD) has been successfully used for many years to provide information about the crystalline structure of samples using coherent scattering of photons from an incident X ray beam. Powder X ray diffraction is commonly used to identify unknown substances by comparing diffraction data against a database [START_REF] Hanawalt | Identification of crystalline materials. classification and use of x-ray diffraction patterns[END_REF][START_REF] Crespy | Energy dispersive x-ray diffraction to identify explosive substances: Spectra analysis procedure optimization[END_REF]. It may also be used to characterize heterogeneous solid mixtures to determine relative abundance of crystalline compounds and, when coupled with lattice refinement techniques such as Rietveld refinement, can provide structural information on unknown materials. Most diffractometers use angular dispersive X ray diffraction (ADXRD) in angular dispersion of a diffracted monochromatic radiation is used, as it provides more precise diffraction spectra(see [START_REF] Langford | Powder diffraction[END_REF] for a review). As detailed in ref. [START_REF] Langford | Powder diffraction[END_REF] many deconvolution and pattern fitting procedures were designed for ADXRD data.

However ADXRD presents some drawbacks as it requires most of the times a rotation of the detector and the sample. When no movement is possible in the setup (e.g. as in security screening) or when crystalline structure is dynamically changing (to monitor the change of lattice parameters with pressure or temperature), one must use Energy dispersive X ray diffraction (EDXRD) for which the scattering angle is fixed. In EDXRD, a collimated polychromatic beam of X-rays falls on the sample, and the spectrum of the radiation scattered at a fixed angle θ is analyzed. The spectrum thus obtained is material specific since it is linked to the atomic planar spacing (d) and to the radiation wavelength (λ) according to Bragg law:

2 d sin θ 2 = n λ , (1)
where n is an integer. Using an ideal diffractometer with perfect collimation and detectors, this diffraction spectrum would be composed only of lines predicted from Bragg theory.

As they provide intense radiation with known spectrum and polarization, most EDXRD experiments use synchrotron radiation [START_REF] Buras | X-Ray Energy Dispersive Diffraction[END_REF][START_REF] Buras | X-ray energydispersive diffractometry using synchrotron radiation[END_REF]. Furthermore as angular aperture of synchrotron light is very small, the blurring of the spectra caused by the width of possible diffraction angles is negligible compared to the blurring effect due to the detector response [START_REF] Glazer | Energy-dispersive powder profile refinement using synchrotron radiation[END_REF]. As a consequence most pattern fitting algorithms [START_REF] Albinati | The rietveld method in neutron and x-ray powder diffraction[END_REF][START_REF] Dong | A full-pattern fitting algorithm for energy-dispersive x-ray diffraction[END_REF][START_REF] Honkimaki | Energy-dispersive diffraction with synchrotron radiation and a germanium detector[END_REF] are designed for EDXRD data with synchrotron radiation and consider a spectraly invariant Gaussian peak-shape function [START_REF] Glazer | Energy-dispersive powder profile refinement using synchrotron radiation[END_REF][START_REF] Albinati | The rietveld method in neutron and x-ray powder diffraction[END_REF].

However if synchrotron is an ideal source for EDXRD, it is scarcely available and only X-ray tube can be used in laboratory or for in situ applications. Angular aperture of collimated X-ray source is quite large and consequently geometric blurring effects result in a "smoothing out" of the perfect line spectra into a less well resolved continuous spectra [START_REF] Luggar | Determination of the geometric blurring of an energy dispersive x-ray diffraction (EDXRD) system and its use in the simulation of experimentally derived diffraction profiles[END_REF]. Very tight collimation, whilst minimizing this effect, would make inefficient use of the scattered flux available and then drastically increase the acquisition time.

As geometric effect introduces an energy variant blur, such EDXRD cannot be processed by existent algorithms developed for ADXRD or EDXRD with synchrotron source. In this paper we propose an inverse problem based method for restoration of EDXRD diffraction spectra and direct extraction of diffraction peaks from crystalline material. It uses a physical model of the measurements de-

Detector

Object

Source

x y z scribed in section 3. Then in section 4, following an inverse problem approach, we introduce some priors and define the maximum a posteriori (MAP) solution which is estimated (in Sec. 5) using a continuous optimization procedure. Finally, we apply our method for estimation of inter-planar spacing of real EDXRD spectrum acquired in laboratory.

EDXRD setup and notation

As depicted in figure 1, the EDXRD setup consists of a collimated poly-chromatic X-ray source and a collimated spectroscopic detector D placed at nominal angle θ of the beam axis. The studied object O is at the center of the diffractometer defined as the intersection of the scatter collimator axis and the source collimator axis.

Let's suppose that we observe at angle θ a pure crystalline powder sample under ideal conditions: using a source with flat spectrum, a diffractometer with an infinite angular resolution, a perfect detector without blur or noise. Without any other spurious effects like Compton scattering, coherently scattered photons are detected only at specific wavelengths related to the atomic interplanar spacing d of the observed material according to the Bragg law. The ideal diffraction spectrum p(ξ) will thus be composed of only few Bragg peaks specific of the material crystalline structure.

Unfortunately, in practice, several effects alter the EDXRD spectrum leading to erroneous detections by fur-ther material identification algorithms. The purpose of the method presented in this paper is to recover the ideal diffraction spectrum p(ξ) from such altered measurements y by means of an inverse problems approach. To this end, we propose an linear modelization of the successive degradations leading to the data y. Using matrix notation ( • is matrix multiplication and × element wise multiplication), this models writes:

y = B • t × (H • p + g) + n , (2)
where:

H is the operator that models the blur caused by the angular distribution of the system.

g is an additional baseline spectrum caused by the amorphous part of the sample incoherent scattered radiation and multiple scattering.

t accounts for the attenuation within the sample, the non flat emission spectrum of the X ray tube and the fluorescence of the anode and the collimators.

B is the blurring operator that models the finite energy resolution of the detector.

n accounts for photon counting and detector noises.

This modelization and corresponding approximations are detailled in the next section.

3. The direct problem: from perfect spectrum to real measurements

Alterations caused by the setup

We suppose that the surface defined by the source focal spot and the detector active region are orthogonal to the z axis at fixed depth of s z and d z respectively. As described on figure 1

θ = acos   -→ so . -→ od -→ so -→ od   . (3)
Owing to the finite angular resolution of the setup, this angle may be different from the nominal angle θ defined by our diffractometer. This is the cause of a blur on the spectrum. Indeed, according to the Bragg's law, a diffraction peak located at an energy ξ for the nominal angle θ will be observed at the energy E for any angle θ according to the following relation:

E = ξ sin( θ 2 ) sin( θ 2 ) . (4)
Conversely, for any angle θ, the measured intensity at the energy E is proportional to the intensity of the nominal diffraction spectrum p E sin(θ/2) sin( θ/2) . As consequence, after integration over the source focal spot surface S, the object volume O and the active region of the detector D, the spectrum f (E) measured by a perfect detector is:

f (E) = D O S N (E, s, o) A(E, s, o, d) C s (s, o) . . . C d (o, d) L(θ) G(s, o, d) p E sin(θ/2) sin( θ/2) ds do dd . ( 5
)
where: • L(θ) is the Lorentz polarization factor [START_REF] Guinier | Theorie et Technique de la Radiocristallographie[END_REF];

• N (E,
C i (a, b) = 0 if (ab) is obstructed by the collimator, 1 else, (6)
• G(s, o, d) is a geometrical factor depending only on distances between the setup elements:

G(s, o, d) = 1 -→ so 2 -→ od , (7)
In some cases, if the object is not too extended, it is possible to assume that Lorentz polarization and geometrical factors are constant. In order to cope with such equation, several simplifications are described in the next subsection.

Attenuation and source spectrum

The source is supposed to emit uniformly over the focal spot and this emission is supposed to be isotropic. Therefore, we can define N (E) the number of photons emitted during the acquisition. When the source is an X ray tube, the source spectrum N (E) is a bremsstrahlung spectrum with fluorescence lines depending of the X ray tube anode.

The attenuation between s, o and d follows the Beer-Lambert law: where µ(E, ℓ) is the linear attenuation of the material at position ℓ at energy E. If the probed volume dimension along the x and y axis are smaller than the morphological structures in the observed object, it can be considered as a succession of layers of different linear attenuation µ i (E) and thickness L i orthogonal to the z axis as depicted on figure 2. Moreover, as the scattering angle θ and the angular aperture of collimators α are sufficiently small, it is possible to hold the classical small angle approximations cos(θ) ≈ 1 and as proposed by Harding [START_REF] Harding | Status and outlook of coherent-xray scatter imaging[END_REF][START_REF] Delfs | Energy-dispersive coherent scatter computed tomography[END_REF][START_REF] Pani | Modelling an energy-dispersive x-ray diffraction system for drug detection[END_REF]. The attenuation term thus simplifies to the transmission attenuation:

A(E, s, o, d) = exp - o s µ(E, ℓ) dℓ - d o µ(E, ℓ) dℓ , (8)
A(E) = exp - i µ i (E) L i . (9)
With such approximations, the product t(E) = N (E) A(E) represents the number of transmitted photons which may be estimated using a spectroscopic detector placed in transmission direction [START_REF] Harding | Status and outlook of coherent-xray scatter imaging[END_REF].

Model of the geometrical blur

As the transmitted spectrum t(E) depends only on the energy, it can be factorized and removed from the integration in equation 5. As consequence, in this integral only the term p E sin(θ/2) sin( θ/2) depends on energy. of the angular distribution of the experimental setup presented at sections 6.1 and 6.3

We define the kernel function h():

h sin( θ/2) sin(θ/2) = Γ(θ) C s (s, o) C d (o, d) L(θ) G(s, o, d) dℓ . (10)
where Γ(θ) is the contour defined by sod = θ and dℓ an elementary part of this contour. This function is a weighted density where h sin( θ/2) sin(θ/2) is length of the iso-contours defined by sod = θ in the seven dimensions volume SOD. These contours are weighted by the geometric, Lorentz polarization and collimators functions.

According to the Bragg law, two couples {E, θ} and {E ′ , θ ′ } probing the same inter-planar spacing d are linked by the equation 4. We thus have:

h sin( θ/2) sin(θ/2) = h E ξ . (11)
The geometric blurring effect caused by the configuration of the setup can be therefore expressed by a simple integral as we can rewrite equation 5:

f (E) = t(E) h E ξ p(ξ)dξ . (12)
We recognize a Fredholm integral equation of the first kind.

Kernel numerical estimation. This equation cannot be analytically solved and we proposed to use a Monte Carlo integration to provide a good numerical estimation of it. Monte Carlo should be less prone to artifacts than Riemann sum proposed by Langford [START_REF] Luggar | Determination of the geometric blurring of an energy dispersive x-ray diffraction (EDXRD) system and its use in the simulation of experimentally derived diffraction profiles[END_REF][START_REF] Pani | Modelling an energy-dispersive x-ray diffraction system for drug detection[END_REF]. This method consists on picking randomly a large number of septuplet

{s x , s y , o x , o y , o z , d x , d y }.
For each septuplet, if both collimator functions are equal to 1, we compute the angle θ and the geometric and Lorentz polarization functions. The density of sin( θ/2) sin(θ/2) weighted by geometrical and Lorentz polarization functions is then estimated using a Parzen-Rosenblatt window method with an Epanechnikov window REFS. An exemple of such histogram computed for experimental conditions given at sections 6.1 and 6.3 is shown on figure 3. This density presents an important negative skewness and cannot be easily approximated by distribution functions classically used in ADXRD peak fitting algorithms (such as Gaussian, Cauchy, Pearson. . . ) [START_REF] Langford | Powder diffraction[END_REF].

Background

During the acquisition, the detector measures not only blurred diffraction peaks but also an additional baseline spectrum. This spurious background g(E) level is caused by amorphous parts of the sample, incoherent scattered radiation and multiple scattering.

m(E) = t(E) (f (E) + g(E)) . (13)

Acquisition

The detector measures the number of photons in a discrete number of energy channels. These measurements y are inevitably corrupted by some noise n. There are mainly three sources of noise:

• photon counting noise which follows a Poisson law,

• thermal noise which is a Gaussian noise,

• modelization noise which is use to take into account approximations made to define the physical direct model. We suppose this noise Gaussian.

For a channel j centered in E j , the detector has a spectral response b (j) (E) and the measurement y j in this channel is:

y j = t(E) (f (E) + g(E)) b (j) (E)dE + n j , (14)
If the response is the same for all channels (b (j

) (E) = b(E -E j ))
, the equation 14 is a convolution. In practice the support of b (j) (E) is larger than the energy step between two consecutives channels, as consequence it causes spectral blurring.

Discretization

Without loss of generality, the spectral distributions t(E), g(E) and p(ξ) can be parametrized by a finite number of parameters by means of expansion onto a basis of interpolation functions (REF). As the measurements y are by definition discrete, it is natural to use the same evenly spaced spectral grid G E = {E ℓ ; ℓ = 1, . . . , N E } to discretize t(E), g(E) and the spectral response functions b (j) (E). In this case, there is no need to explicitly choose any particular basis of interpolation functions as far as it has the properties required to make them suitable as interpolation functions and it comes:

t i = t(E i ) (15)
g i = g(E i ) (16) b (j) i = b (j) (E i ) (17)
For the parameter p(ξ) which is the one we want to extract, more attention has to be paid on discretization. It can be discretized on another spectral grid G ǫ = {ǫ j ; j = 1, . . . , N ǫ } using the interpolation function v j (ǫ)

p(ξ) = j p j v j (ξ-ǫ j ) . ( 18
)
The grid G ǫ can be the same as G E but can also be more finely discretized to perform super-resolution. Furthermore, it may not be necessary evenly spaced and can be designed to be evenly spaced in the momentum transfer space (

x j =ǫ j sin( θ/2) h c
) or in the angle space for a given energy Ē in order to be directly compared with a existing database indexed in momentum transfer or in angle. In this work, we choose to use the grid of the measurements:

G ǫ = G E (ǫ j =E j ) and to define v j (ǫ) = v(ǫ) as uniform: v(ǫ) = 1. κ if -κ 2 ≤ ǫ < κ 2 0 otherwise, ( 19
)
where κ is the channel width of the detector.

Using such parametrization our discrete model writes:

y i = j B i,j t j g j + k H j,k p k + n i . ( 20
)
where the coefficients of B i,j and H j,k are:

B i,j = b (i) (E j ) (21)
H j,k = R + * h Ej E ′ v(E ′ -ǫ k ) dE ′ . ( 22
)
Equation 20 can be written using matrix notation and becomes the direct discrete:

y = B • t × (H • p + g) + n . ( 23
)
This equation defines our direct problem and thus can be used in order to simulate the data. In practice filter reponse b (i) remains the same for all channels and the convolution matrix B can be fastly applied by means of a Fourier transform. A simulation using this model with g = 0 and p = 0 excepted at 30 keV, 80 keV and 140 keV is shown on figure 4 for experimental conditions given at sections 6.1 and 6.3. On this figure, it can be clearly seen that the spread of the peaks increases as their energy rises.

The inverse problem:

The goal of the presented work is to estimate the diffraction peak distribution p from the data y according to the model defined by equation 2. We suppose that the operators B and H are known and that a good estimation of t is provided by measurements from a spectroscopic detector placed in transmission. However, even if it is not of interest, we still have to estimate the background g jointly with p (both of size N ) from a data spectrum of size N . This is a typical inverse problem, that can be solved in a penalized likelihood or maximum a posteriori (MAP) framework [START_REF] Thiébaut | Strict a priori constraints for maximum likelihood blind deconvolution[END_REF]. This is achieved by estimating the couple {p + , g + } that minimizes the cost function ξ(f , g):

{p + , g + } = arg min {p,g} ξ(p, g), (24)
ξ(p, g) = Φ(p, g, y) + Ψ(p) + Θ(g), (25)
This cost function ξ(p, g) is the sum of three terms: a likelihood penalty Φ(p, g; y) ensuring the agreement between the model defined by equation 2 and the data y, and two regularization penalties Ψ(p) and Θ(g) introducing subjective a priori knowledge about the diffraction peak distribution p and the background g respectively.

The likelihood penalty term

Assuming Gaussian noise, the likelihood penalty reads:

Φ(p, g) = r T • W • r , (26)
with the residuals:

r = y -B • t × (H • p + g) , (27)
and where the weighting matrix W is the inverse of the covariance of the noise: W = C -1 noise . If the noise is uncorrelated, the weighting matrix is diagonal W = diag(w) and Eq. 26 simplifies to:

Φ(p, g) = i w i r 2 i , (28)
where 1/w i is the noise variance in the channel i. This model can cope with non-stationary noise and can be used to express confidence on measurements on each spectral channel of the data. Since unmeasured data have infinite variance, we can readily deal with missing (outside the field of view) or bad channels as follows:

w i ≡ Var(y i ) -1 if y i is measured, 0 otherwise. ( 29
)
Except for very low detector noise, we can approximate the total noise (Gaussian detector noise plus Poissonnian signal noise) by a non stationary Gaussian noise [START_REF] Mugnier | Mistral: a myopic edgepreserving image restoration method, with application to astronomical adaptive-optics-corrected long-exposure images[END_REF]:

w i ≡ γ max(y i , 0) + σ 2 i -1 if y i is measured, 0 otherwise, (30)
where γ accounts for the quantization factor of the detector and σ 2 i is the variance of other approximately Gaussian noise (for example read-out noise) in the channel i.

Priors

The most evident objective prior on interference functions we can introduce is its positivity. Unfortunately this constraint is not strong enough to lead to acceptable solutions, and other a priori knowledge are needed. These priors are introduced by the use of the regularization functions Ψ(p) and Θ(g). As the diffraction peak distribution p and the background g have very different a priori morphological shape, these two functions are designed separately.

Prior on the diffraction peaks distribution. Interference function of a crystalline material is non-null only for values that correspond to the inverse of crystal inter-planar spacing. The diffraction peak distribution p is then very sparse. Recovery of sparse vectors has been extensively studied in this last decade 1 . Donoho [START_REF] Donoho | For most large underdetermined systems of equations, the minimal l1-norm near-solution approximates the sparsest near-solution[END_REF] has shown that minimizing the ℓ 1 norm of a vector p is the best convex regularization function that favor sparsest solutions. Unfortunately, ℓ 1 norm is not continuously derivable preventing fast convergence of our continuous optimization procedure. For that reason we use a continuously derivable hyperbolic approximation defined by:

Ψ(p) = k α k p 2 k + ε 2 -ε , (31)
where ε is very small (≈ 10 -9 ) and α k are the so-called hyper-parameters that may vary with the channel k.

Priors on the background. As we don't have precise (i.e. quantitative) information about the main effects that lead its formation, we propose to introduce some smoothness priors about the background g. The background regularization function Θ(g) is thus a classical quadratic smoothness regularization function defined by:

Θ(g) = k β k (g k+1 -g k ) 2 , (32)
where β k are channel varying hyper-parameters. This regularization function will act as a low pass filter which prevent sharp changes of the interference function.

Hyper-parameters settings

α k and β k are the hyper-parameters that are used to balance the influence of both priors and likelihood and that have to be properly tuned. Unfortunately, as these hyper-parameters are channel varying, we replace an estimation problem with 2 N parameters by another problem of (2 N -1) hyper-parameters estimation. However, the large number of hyper-parameters defined in this very general scheme can be drasticly reduced according to simple rules.

In this work, the parameter of interest is the diffraction peak distribution p. Without informations about the material nature, the probability of presence of a diffraction peak is uniformly distributed over the channels. Unfortunately, constant hyper-parameters α k and β k favor diffraction peaks at low energy. Indeed, as shown in figure 4, in the data y, diffraction peaks at low energy are sharper than those at high energy and a sharp peak has a higher probability to be correlated with some noise realization than a smooth one. The hyper-parameters α k have thus to be tuned to enforce an uniform distribution of the diffraction peaks. Mimicking the L2 normalisation of atoms of a dictionnary used in matching pursuit algorithms [START_REF] Mallat | Matching pursuits with time-frequency dictionaries[END_REF], we propose a normalization of the hyper-parameter α k by the L2 norm of the corresponding blurred peaks. It writes thus:

α k = α i   j B i,j H j,k   2 , ( 33
)
where α is a global hyper-parameter for diffraction peak distribution regularization. Furthermore, the background regularization Θ(g) favors smooth background. As a peak at high energy is smooth, the hyper-parameters β k must be tuned to enforce a background smoother than diffraction peaks everywhere in the spectrum and to enforce an uniform distribution of the distribution peaks.

β k = β 1 i [∇ ( H -,i )] 2 k , ( 34
)
where H -,i means the i-th column of H and ∇ is the finite difference operator. β is a global hyper-parameter for the background regularization. This normalization of β k ensure that the regularization function Θ(g) gives the same penalty for all peaks wherever their position is in the spectrum. Finaly only two hyper-parameters β and α have to be tuned for the reconstruction.

Algorithm summary

In our framework, the reconstruction of the diffraction peak distribution p and the background g corresponds to the determination of the couple {p + , g + } that minimizes the objective criterion defined in Eq. 25. The optimal reconstructed diffraction peak distribution depends on the particular values of the hyper-parameters α and β. Choosing their optimal values is difficult and it is done for now by trials and errors. However, for future automatic hyperparameters estimation, the L-curve method [START_REF] Hansen | Analysis of discrete ill-posed problems by means of the l-curve[END_REF] or Stein unbiased risk estimator [START_REF] Stein | Estimation of the mean of a multivariate normal distribution[END_REF] seems suitable, but this deserves an extensive study which is out of the scope of this paper.

To minimize the criterion defined in Eq. 25., we use the VMLM-B (Variable Metric with Limited Memory-Bounded) algorithm [START_REF] Thiébaut | Optimization issues in blind deconvolution algorithms[END_REF] which is a limited memory variant of the variable metric method with BFGS updates [START_REF] Nocedal | Numerical optimization[END_REF]. This algorithm can further account for bound constraints on the parameters. We make use of this feature to enforce positivity of p and g. This algorithm has proved its effectiveness for image reconstruction and only requires the computation of the penalty function being minimized and its gradient. The memory requirement is a few times the size of the problem.

Results

In this section, we illustrate the performance of the proposed algorithm using real data. Two crystalline materials with well known inter-planar spacing (since [START_REF] Bragg | The structure of some crystals as indicated by their diffraction of x-rays[END_REF][START_REF] Hull | A new method of x-ray crystal analysis[END_REF]) are used in these experiments: NaCl and graphite.

Experimental setup

Schematic representation of the experimental setup is shown in figure 5. The sample is illuminated by a 140 kV tungsten anode X-ray tube operated at 280 µA and placed at 36.7 cm from the sample. Its focal spot is 1 mm wide. The incident X-ray beam is collimated by a tungsten slit of 0.2 mm in width, 8 mm in height and 8 cm in length. This source collimator is placed 18.4 cm away from the focal spot.

The scattered photons are detected with an high purity germanium (HPGe) detector oriented at θ = 4.26 • and placed 24 cm from the sample. The detector collimator located in front of the detector is a slit of 8 mm in height and 8 cm in length. This slit has a width of 0.02 mm on the detector side and 0.2 mm on the other side. The spectral response function of the HPGe detector energy b is supposed to be a Gaussian with a full width at half maximum (FWHM) of 0.4 keV. It is supposed to be constant on the whole energy range considered in these experiments: b (j) = b. We consider every photon detected in N L = 1818 energy channels of the detector between 15 keV and 160 keV (sampling step size of 80 eV). With such a couple source/detector the range of the probed interplanar spacing is beetwen 0.2 nm -1 and 2.4 nm -1 . The exposure time is set at 600 s, for each spectra to have a good counting statistic.

In both experiments, the matrix H is computed according to the equation 22 using these setup parameters and the approximative shape of the sample. Without transmission detector, the source spectrum was approximated by a Birch and Marshall [START_REF] Birch | Computation of bremsstrahlung x-ray spectra and comparison with spectra measured with a ge (li) detector[END_REF] semi empirical spectrum for a tungsten anode 160 kV source. For the attenuation spectrum we use the theoretical attenuation of the samples.

NaCl

The NaCl used was held in a thin walled polyethylene cylinder of 1 cm diameter and 2 cm in height. The scattering measurements are presented on figure 6. Several comments can be made on these data and we can notice:

• a continuous background that is higher for low energy,

• four peaks that may be produced by diffraction at 50 keV, 59 keV, 83 keV and 102 keV,

• two peaks at 59.7 keV and 67 keV that are caused by tungsten fluorescence,

• a very broad peak between 30 keV and 40 keV which may be due to the scattering of the plastic cylinder. This peak is too broad to be caused by scattering of crystalline materials. These data were processed using the proposed algorithm with hyper parameters α = 0.05 and β = 10. The source spectrum was supposed to follow the Birch and Marshall synthetic spectrum and we use the theoretical attenuation of a 1 cm thick NaCl cube. Estimated diffraction peak distribution p and background g are shown on figure 7.

As expected the estimated background g fits the broad peak between 30 keV and 40 keV but not the next peak as 50 keV which is a diffraction peak. As the background plotted on figure 7 is the product g × t, the peaks about 58 keV and 67 keV correspond to the Tungsten fluorescence lines.

The diffraction peak distribution p shows six peaks. The positions of the centroid of these peaks are presented in the table 1. They can be compared with the theoretical diffraction peak positions estimated for θ = 4.26 7) compared to position for θ = 4.26 • estimated using inter-plannar spacing measurements from from [START_REF] Swanson | Standard x-ray diffraction powder patterns[END_REF]. false is for false detection and ND for non detection. The theoretical position of the peaks have been estimated for θ = 4.26

•
ASTM card 5-0628). Every peaks present a very accurate position compared to the Swanson estimation with an error smaller than 1 keV. The last detected peak (at 131.29 keV) presents the higher positioning error which can be explained by the relatively small number of photons measured around this energy. The five diffraction peaks that were not detected, are the five peaks with the lower relative intensity (< 10%) in the measurements of Swanson.

A small peak has been falsely detected at 64.1 keV. This false detection may be due to the tungsten fluorescence line K β at 67.25 keV. Indeed, as we use only a very coarse estimation of both source spectrum and material attenuation, the magnitude of these rays may be under-estimated, and their residuals may be explained as a scattering line by our algorithm. This problem of too coarse estimation of source spectrum and material attenuation is responsible for the misestimation of peak magnitude. Furthermore as the presented method is for powder diffraction, multiplicity and texture (prefered orientation of crystal) may affect measured relative intensity. As consequence, only peak's positions can be used for material identification and not relative intensity.

The figure 8 shows the model of the data computed with estimated solutions p and g. This model shows a good accordance with the data estiblishing the validity of our linear model of the formation of EDXRD data. However, it can be noticed that the model is always lower than data around the diffraction peaks. This under estimation can be explained by a too coarse estimation of the setup geometry that has lead to an over estimation of the angular dispersion giving thus too broad diffraction peaks.

Graphite

The graphite used in this experiment was a 4 cm thick piece of about 2 cm in width and height. This large size of the sample is responsible for quite broad scattering peaks in the scattering measurements. This EDXRD diffraction spectrum is presented on figure 9. Data were processed using the proposed algorithm with hyper parameters α = 4.10 -3 and β = 5. We use the theoretical attenuation of a 4 cm thick graphite piece. The results of the algorithm is plotted on figure 10 and the model of the data computed with estimated solutions p and g is plotted on figure 11. This figure shows the accordance between the model and the data. The estimated locations of diffraction peaks are shown in the table 2 and can be compared with "theoretical" peaks position estimated using inter-planar spacing measurements given by ASTM card 12-212 [START_REF] Read | ASTM card 12-212, bell Telephone Laboratories[END_REF].

On figure 10, we can see that seven diffraction peaks are detected. Most of them (excepted the highest peak at 49.11 keV) were difficult to distinguish on the data (figure 9), in particular for those at energy higher than 120 keV where very few photons (< 5) are detected in a single channel. Only the three fainter peaks were not detected. Furthermore, the K α line of the tungsten (at 58.9 keV) is not identified as a diffraction peak. The highest peak at 49.11 keV (hkℓ : 002) has two visible harmonics at 98.50 keV (hkℓ : 004) and 149.47 keV (hkℓ : 006). The baseline of the diffraction spectrum is well fitted by the smooth background term g.

Discussion

In this study we described a numerical method to restore energy dispersive X-ray diffraction data and to extract diffraction peaks. It is based on an inverse problem approach. It uses a linear direct model derived from an analysis of the formation of the data and the causes of diffraction peaks blurring. Using a continuous optimization technique, it then estimates the "perfect" diffraction spectrum that best reproduces the data according to the model and some additional priors. These priors are the sparsity of the crystalline diffraction spectrum and the smoothness of the baseline spectrum.

The level of influence of the priors are tuned using two hyper-parameters. In the present work the hyperparameters were set by trials and error for each spectrum. The setting of hyperparameters can be tedious and consequently it will be necessary to find an automatic method for routine use of this method. A further study will be necessary to determine whether existing hyper-parameters setting methods (L-curve [START_REF] Hansen | Analysis of discrete ill-posed problems by means of the l-curve[END_REF], GCV [START_REF] Golub | Generalized cross-validation as a method for choosing a good ridge parameter[END_REF], SURE [START_REF] Stein | Estimation of the mean of a multivariate normal distribution[END_REF]. . . ) can be used or whether a specific method has to be defined.

This method was tested on real diffraction spectra of known crystalline materials: NaCl and graphite. Computation time were approximately equal to acquisition time (several minutes). On both materials, the presented method was able to extract precisely diffraction peaks even at high energy where the number of detected photons per channel is low (<10). These real results clearly assess two facts: (i) the direct model is correct as it provides good fits of the data and (ii) the priors effectively provide satisfying solutions.

As attenuation and source spectra were only coarsely estimated, the estimated relative magnitude of the diffrac-tion peaks is uninformative for material identification and only peaks location could be used. This coarse estimation is also responsible of false detections around fluorescence line of the tungsten. These drawback should be attenuated by the simultaneous acquisition of the transmission spectrum that should give (in the small angle approximation) a better estimation of the attenuation times source spectrum.

The other sticking point of the presented method is the necessity to have at least a coarse estimation of the shape of the sample to compute the blurring kernel h(). A way to bypass this shape estimation step is to directly calibrate this kernel h() using the measurements as it is done in image blind deconvolution algorithm [START_REF] Soulez | Blind deconvolution of video sequences[END_REF].

A natural extension of this "inverse problem" is the reconstruction of a 3D volume using an array of detectors and the estimation of a whole 3D crystalline material as we proposed in [START_REF] Soulez | Spatial super-resolution for line parallel imaging in energy dispersive x-ray diffraction[END_REF].

Figure 1 :

 1 Figure 1: Schematic view of a X ray diffraction setup and notations used in the model.

  consider a triplet {s, o, d} composed of three points s = [s x , s y , s z ] T on the source focal spot, o = [o x , o y , o z ] T in the object and d = [d x , d y , d z ] T on the active surface of the detector. For this triplet, the scattering angle θ is defined by:

  with i = s and (a, b) = (s, o) or i = d and (a, b) = (o, d)

1 Figure 2 :

 12 Figure 2: Schematic view of a X ray diffraction setup and notations used in the model.

Figure 3 :

 3 Figure 3: Density h sin( θ/2) sin(θ/2)

Figure 4 :

 4 Figure 4: Total response of the setup for three peaks of equal intensity placed at 30 keV, 80 keV and 140 keV.

Figure 5 :

 5 Figure 5: Schematic view of the setup used in the experiments.

Figure 6 :

 6 Figure 6: Measurements on a NaCl sample.

Figure 7 :

 7 Figure 7: NaCl diffraction spectrum (in black) with the coresponding estimated background g × t (in dashed light gray) and diffraction peaks p (in dark gray).

Figure 9 :

 9 Figure 9: EDXRD measurements on a Graphite sample.

Figure 10 :Figure 11 :

 1011 Figure 10: Graphite diffraction spectrum (in black) with the coresponding estimated background g × t (in dashed light gray) and diffraction peaks p (in dark gray).

Table 1 :

 1 Estimated diffraction peak position (shown on Fig.

	• us-

Table 2 :

 2 Estimated and theoretical position of graphite diffraction peaks shown on figure 10, in keV. * denotes two harmonics of the peak at 49.11 keV. false is for false detection and ND for non detection.

					26 • peak position
			in %	in keV	in keV
	002	3.370	100	49.49	49.11
	100	2.132	2	78.23	77.66
	101	2.036	3	81.92	80.69
	102	1.800	1	92.66	ND
	004 *	1.682	8	99.16	98.50
	103	1.541	2	108.24	104.76
	110	1.232	6	135.34	133.58
	112	1.155	2	144.41	ND
	006 * -105 1.120	2	148.92	149.47
	201	1.050	2	158.85	ND

The website http://www.compressedsensing.com/ contains many references about sparsity and compressive sensing.

http://yorick.sourceforge.net

Acknowledgment

The algorithms presented in this article have been implemented in Yorick 2 , a free data processing language written by D. Munro. The authors are grateful to Éric Thiébaut whose fruitful advices and attentive reading have helped to improve this paper. This work has been supported by project SPIDERS funded by the French National Research Agency.


Paolo Simonetto

Daniel Archambault

David Auber

Romain Bourqui

ImPrEd: An Improved Force-Directed Algorithm that Prevents Nodes from Crossing Edges

Keywords: Descriptors (according to ACM CCS): G, 2, 2 [Discrete Mathematics]: Graph Theory-Graph Algorithms

is a force-directed algorithm that improves the existing layout of a graph while preserving its edge crossing properties. The algorithm has a number of applications including: improving the layouts of planar graph drawing algorithms, interacting with a graph layout, and drawing Euler-like diagrams. The algorithm ensures that nodes do not cross edges during its execution. However, PrEd can be computationally expensive and overlyrestrictive in terms of node movement. In this paper, we introduce ImPrEd: an improved version of PrEd that overcomes some of its limitations and widens its range of applicability. ImPrEd also adds features such as flexible or crossable edges, allowing for greater control over the output. Flexible edges, in particular, can improve the distribution of graph elements and the angular resolution of the input graph. They can also be used to generate Euler diagrams with smooth boundaries. As flexible edges increase data set size, we experience an execution/drawing quality trade off. However, when flexible edges are not used, ImPrEd proves to be consistently faster than PrEd.

Introduction

PrEd [START_REF] Bertault | A force-directed algorithm that preserves edge crossing properties[END_REF] is a force-directed algorithm [START_REF] Eades | A heuristic for graph drawing[END_REF][START_REF] Fruchterman | Graph drawing by force-directed placement[END_REF] that preserves the edge crossing properties of an input layout. The algorithm computes the maximal displacement of nodes and restricts their movement to guarantee this property.

The author of PrEd suggested two possible applications: improving the layout [START_REF] Purchase | Which aesthetic has the greatest effect on human understanding?[END_REF] of planar, straight-line graph drawing algorithms [START_REF] De | How to draw a planar graph on a grid[END_REF] and driving force-directed graph layout interactively. By realising that not only does PrEd preserve all edge crossings in an existing layout but actually prevents nodes from crossing edges, we can apply this algorithm to a wider range of problems beyond those stated in the original paper. For example, PrEd can be used to improve Euler-like diagrams [START_REF] Simonetto | Fully automatic visualisation of overlapping sets[END_REF] where nodes are set elements and edges are set boundaries. In general, we can use PrEd for tasks where we would like to bound portions of a layout, extending the applicability beyond "preserving edge crossing properties" [START_REF] Bertault | A force-directed algorithm that preserves edge crossing properties[END_REF].

The contribution of this paper is ImPrEd: an improved force-directed algorithm that prevents nodes from crossing edges. ImPrEd is faster than PrEd and grants greater output control. Our test cases provide some evidence that it has improved visual quality. The modifications made, initially motivated by the problem of drawing Euler-like diagrams, can also be useful for general graph drawing.

Previous and Related Work

Many approaches aim to constrain or restrict node movement in force-directed algorithms. Mental map preservation techniques, described in Section 2.1, restrict node movement for dynamic graph drawing. Section 2.2 describes constraintbased techniques, which offer guarantees on properties in the output of the graph. Section 2.4 discusses planar graph drawing -an application area of PrEd. Finally, PrEd is described in Section 2.5.

Mental Map Preservation

In the field of dynamic graph drawing, a number of forcedirected algorithms have included in their objective functions ways to keep similar nodes in similar areas of the plane. submitted to Eurographics / IEEE Symposium on Visualization 2011 (EuroVis 2011)

Erten et al. [EHK * 04] considered inter-timeslice edges that link the same node across timeslices to drive the positions of nodes to certain area of the plane. Weighting functions have been used to drive, but not constrain, the movement of nodes [START_REF] Frishman | Dynamic drawing of clustered graphs[END_REF][START_REF] Frishman | Online dynamic graph drawing[END_REF]. The PIGALE [START_REF] De | Public implementation of a graph algorithm library and editor (PI-GALE)[END_REF] library also includes a spring embedder which preserves the planar map. Additionally, certain force-directed methods have impeded the movement of nodes via geometric restriction [START_REF] Saffrey | The "mental map" versus "static aesthetic" compromise in dynamic graphs: A user study[END_REF]. For energy-based approaches, penalising terms can be added to discourage nodes from moving too far from their previous positions [START_REF] Boitmanis | Visualizing internet evolution on the autonomous systems level[END_REF].

These approaches modify their respective algorithms to drive nodes in the layout towards certain positions. However, they do not guarantee that a node will stay in a particular region of the plane or that it will not cross edges.

Constraint-Based Graph Layout

Several approaches use hard constraints to guarantee that certain properties of an initial graph layout are not broken [START_REF] Dwyer | Dig-CoLa: Directed graph layout through constrained energy minimization[END_REF][START_REF] Dwyer | IPSep-CoLa: An incremental procedure for separation constraint layout of graphs[END_REF][START_REF] Dwyer | Topology preserving constrained graph layout[END_REF][START_REF] Dwyer | Scalable, versatile and simple constrained graph layout[END_REF]. Properties such as the horizontal or vertical alignment of nodes, nonoverlapping rectangular boundaries, containment of nodes within a boundary, and downward pointing edges in a directed graph can be preserved. The general technique has been shown to be applicable in diverse settings, including edge crossing properties. However, with ImPrEd, we aim at the specific application of node containment which may be handled more simply with a targeted algorithm.

Topology-Driven, Force-Directed Algorithms

Didimo et al. [START_REF] Didimo | Topologydriven force-directed algorithms[END_REF] proposed a framework, combining force-directed and planarisation-based approaches for graph drawing. This hybrid algorithm aims at drawing graphs with a small number of crossings and high angular resolution. Although the algorithm shares similarities with PrEd and ImPrEd, there are significant differences between the approaches. For example, Didimo et al. has forces for geodesic edge tendency which ImPrEd does not. However, this algorithm does not preserve all crossings in the initial layout.

Planar Graph Drawing and Visualisation

Graph visualisation systems are often less interested in strict graph planarity [START_REF] Herman | Graph visualization and navigation in information visualization: A survey[END_REF], since a large graph rarely admits a strict, planar embedding. However, planar graphs are often created to provide an overview of the data. Also, they have recently become of interest to the problem of drawing Euler-like diagrams [START_REF] Simonetto | Fully automatic visualisation of overlapping sets[END_REF]. Current planar graph drawing algorithms [START_REF] De | How to draw a planar graph on a grid[END_REF][START_REF] Gutwenger | Planar polyline drawings with good angular resolution[END_REF] can produce sub-optimal results in terms of graph aesthetics, when considering angular resolution. PrEd has been proposed to improve the output produced by planar graph drawing algorithms, but it can have slow execution time and may produce sub-optimal results for some planar graphs.

Algorithm 1 Pseudo code for the original PrEd approach.

∀v, calculate all F r v , F a v , and F e v ∀v, compute Mv (Figure 1) ∀v, move nodes using min of forces or Mv

An Introduction to PrEd

Like most force-directed approaches, PrEd iteratively refines the layout. Each iteration has three stages. Algorithm 1 gives an overview of the approach. First, a force Fv is computed for every node v ∈ V : using node-node repulsion F r v , edge attraction F a v , and non-incident edge repulsion F e v . Secondly, the algorithm computes the maximal movement permitted Mv for each node as shown in Figure 1. This movement is computed by dividing the directions around each node v into eight sectors and computing the shortest distance between v and the first edge that it could potentially cross. Finally, the positions of the nodes are updated according to the resultant force, ensuring the displacement is less than Mv in the appropriate direction.

PrEd has three forces: node-node repulsion, edge attraction, and node-edge repulsion. The first two forces are:

F r v (u, v) ← δ ||pu -pv|| 2 (pu -pv) F a v (u, v) ← ||pu -pv|| δ 1 (pv -pu)
while the node-edge repulsive force is:

F e v (v, (a, b)) ← (γ -||pv -ve||) 2 ||pv -ve|| (pv -ve) ve ∈ (a, b), a = v, b = v, ||pv -ve|| < γ
In these formulae, px is the position of node x, δ is the optimal distance between two nodes, and γ is the optimal distance between a node and an edge. The value of ve is the projection of node v onto the line defined by the edge (a, b). These forces are computed for every pair of nodes, every edge, and every node-edge pair respectively.

To compute Mv, PrEd divides the 2π angle around each node v into eight congruent angles Av = (A 0 v . . . A 7 v ). The maximal movement is an array of eight values M 0 v . . . M 7 v corresponding to the maximal distance a node can move in the direction of any vector in the sector. Figure 1 depicts these sectors and how node movement is restricted.

Definitions

Before describing the algorithm, we review a few definitions and introduce some notation. A more detailed explanation can be found in graph theory and drawing textbooks [START_REF] Nishizeki | Planar Graph Drawing[END_REF]. where the nodes v ∈ V are assigned to coordinates pv in the plane. The edges e = (u, v) ∈ E are lines connecting nodes at pu and pv and two edges can only cross at their extremities. Some authors refer to this concept as planar graph embedding, whereas a general embedding omits the non-crossing condition. A graph that has an embedding is called planar graph, and a particular embedding is a plane graph.

Let G = (V,
Let G be a plane graph. In its embedding, the edges of the graph divide the plane into regions called faces. Let Φ be the set of faces of G and F ∈ Φ. The boundary B(F) of a face is the set of nodes B n (F) and the set of edges B e (F) that enclose F. The boundary of a face is a collection of nodes and edges for a disconnected graph, a closed walk for a connected graph, and a cycle when the graph is 2-connected. We will use the notation F when F has a connected boundary.

Two embeddings of G are equivalent if the boundary of each face in an embedding corresponds to the boundary of a face in the other. The planar maps of G are equivalence classes of the planar embeddings of G.

A combinatorial embedding of G is a data structure that lists in order the edges incident to the graph nodes. When G is connected, its planar maps correspond to its combinatorial embeddings. When G is not connected, we can obtain nonequivalent embeddings by swapping a connected component of G into a different face, and so the same combinatorial embedding might correspond to different planar maps.

Improvements to the PrEd Approach

ImPrEd is designed to improve the performance and extend the features of PrEd. This section presents the execution time improvements that affect the three steps of Algorithm 1. The surrounding edge computation reduces the number of edges checked when computing the maximal movement of a node. QuadTrees efficiently determine the nodes and edges in a region, reducing the time required for force and maximal displacement computations. New rules for maximal displacement allow nodes to move more freely, enabling faster convergence. Finally, the force system adapts, based on the iteration, leading to improved aesthetics and reliability of input parameters.

Surrounding Edges Computation

The author of PrEd observed that, in certain cases, edges can be crossed only after other edges are crossed. Thus, considering all edges in the maximal movement computation might be redundant. Consider a node v inside a face F of a connected plane graph G = (V, E). By construction, the node must cross an edge of B e ( F) before it can cross any other edge in G. Therefore, we can safely test only the edges of B e ( F). Moreover, these edges are the most influential in the node-edge repulsion force F e v and can be the only edges considered without significantly affecting the final layout. These edges, the surrounding edges (Sv) of the node v, will be the only edges considered in the calculation of F e v and Mv in Algorithm 1. Since the algorithm preserves layout topology, Sv can be computed in pre-processing step and does not need to be updated afterward.

Computing Sv requires us to analyse the plane graph as shown in Figure 2. For now, assume that the graph G is a plane graph, but we will later extend this result to non-plane graphs (a non-plane graph can either be a non-planar graph or a planar graph with a non-planar embedding). The surrounding edges of a node v ∈ V can be computed as the union of all edges in the faces that contain v (see Figure 2c):

Sv = F∈Φv B e (F) where Φv = F ∈ Φ | v ∈ B n (F)
In a general plane graph, the boundary of F might not be connected and can depend on the embedding of the graph. We compute the connected components and the boundary of their faces F. As we are dealing with connected components, we can compute their faces using their combinatorial embedding. Finally, we merge them to obtain the boundary of the face F (see Figure 2d).

We decompose the graph into its connected components and compute the faces expressed by the embedding. Each component, C i , has a single unbounded or external face, F ext i , and a set F int i of zero or more bounded or internal faces F int i, j . Moreover, as there are no edge crossings, each connected component is contained in exactly one face.

We construct a partial ordering between components, so that the components C i and C j are in the relationship C i < C j if and only if C i is contained by an internal face of C j and C j is contained by the external face of C i . The partial ordering organises the components into a hierarchy. The resulting forest has unrelated components as roots and all other components ordered by containment (see Figure 2b). Let us assume that G is not a plane graph. To apply the above-described method, G is reduced to a plane graph, Ĝ, through planar augmentation [START_REF] Battista | Graph Drawing; Algorithms for the Visualization of Graphs[END_REF]. We label all nodes and edges of Ĝ using elements that generated them, and these labels are used to compute Sv. The nodes of Ĝ that do not correspond to a node in G are not labeled and are not considered as part of the face boundary (Figure 3).

The complexity of the surrounding edges computation depends on the initial layout configuration. If G is a plane graph, each edge is considered twice. Since the number of edges and faces is linear with respect to the number of nodes, the overall complexity is O(|V |). If G is not a plane graph, in worst case, we have O(|E| 2 ) crossings. As each crossing produces at most one node and at most two segments and outputs a plane graph, the overall worst case complexity is O(|V | + |E| 2 ). In general, the average number of surrounding edges is much smaller |E|. However, there are cases (for example, trees) where this worst case is realized. It is im-portant to note, however, that this step is pre-processing and does not influence the main cycle of the algorithm.

QuadTrees

We use QuadTrees to determine with sublinear complexity the nodes N d v and edges E d v at a distance d from a node v. The data structure has been used in force-directed algorithms [START_REF] Quigley | FADE: Graph drawing, clustering, and visual abstraction[END_REF] as distant elements can often be approximated or ignored. In ImPrEd, when computing F r and F e , we only consider nodes in N 3δ v and the edges

E γ v .
The set of nearby edges can also be used to reduce the running time of the maximal movement computation. For stability reasons, ImPrEd defines a global maximal distance d that no node can exceed when moving (see Section 4.4). Consider a node v and an edge e that cannot be crossed. If v and the extremities of e move closer to each other, they cross only if v and e are closer than 2 d. We can then safely consider only the edges Sv ∩ E 2 d v in the computation of Mv.

The edge sets complement each other, overcoming many of their individual limitations. For example, when using Im-PrEd on a tree, we have Sv = E as all edges are surrounding edges for every node. In this case, E 2 d v helps alleviate the problem. Conversely, when using ImPrEd on a graph with an unbalanced layout, Sv often keeps the number of tested edges small.

New Rules for the Maximal Movement Calculation

PrEd guarantees that no new crossings can be incurred or undone. To enforce this property, the algorithm considers each node-edge pair (v, e), v ∈ V , e = {w, z} to ensure that v, w and z do not create new or undo crossings. which is perfectly safe. Moreover, a node inside an acute angle or between two edges will have its movement unnecessarily constrained. By allowing nodes to escape these highstress configurations, we are likely to obtain a stable configuration quickly, with additional positive effects on drawing quality. Therefore, we modify the node restriction rules to maximise node movement, while preserving the crossings in the initial layout. Figure 4 and 5 compare the movement sectors computed with the old and the new rules. This improvement affects the computation of Mv in Algorithm 1.

Each node v has an associated collision vector cv. The vector is the minimum movement of v before crossings are incurred. Let ve be the projection of v in the line where e lies. The collision vector cv is defined as:

Case ve ∈ e : cv = - vve Case / ∈ : cv = → vw, -vz)
The node can most cv|| in without crossing e. both and z can move as well. sider a line l perpendicular to and crossing its mid-We redefine cv and define cw and cz as the vectors connecting w and z to their projections on l. Again, if nodes move less than their collision vectors, they will not cross l and therefore no crossing between e and v can occur.

To permit maximal movement, each sector can be extended until it touches, but does not become collinear with, the line l. This result can be achieved by multiplying the length of the collision vector by σ, defined as the secant of the angle between c and closest radius of the sector:

M j x < cx • σ( cx, A j x ) σ( cx, A j x ) =          1 i -j ≡ 0 mod 8 sec cx -( j + 1)π/4 i -j ≡ 1, 2 mod 8 sec cx -jπ/4 i -j ≡ 6, 7 mod 8 ∞ i -j ≡ 3, 4, 5 mod 8
where cx is the angle of the given vector and where i is the index into A i x in which cx lies. In the first case, the vector lies on the sector and the closest radius coincides with cx, thus sec(0) = 1. The second case considers the two sectors in clockwise order from A i

x , where the closest radius corresponds to the sector border of A j x that has angle ( j + 1)π/4. The third case handles the two sectors in anticlockwise order that have the closest border with angle jπ/4. The forth case handles the sectors not facing l that are unrestricted.

To prove the correctness of the algorithm, we first prove that when a node moves it cannot cross an edge. Given a node v and a non-incident edge e, the line l (see Figure 4 and5) divides the plane into two half planes: the first, Rv(v, e), contains v and the second, Re(v, e), contains e = (w, z). For all other λ ∈ E and λ non-incident to v, we compute the half planes that contain v using λ and v . Similarly, for all other u ∈ V and u non-incident to e, we compute the half planes that contain e using e and u. The regions defined by the intersection of these half planes are Rv and Re respectively and are computed for all nodes and edges in the graph. The above-described formula bounds the sectors around every v to be in Rv and every w and z to be in Re because the nodes must remain on the same side of all defining half planes of the region. Therefore, Rv contains any future position of any node v and Re, being convex, as it is the intersection of convex regions [DBCVKO00], fully contains any final position of the edge e. Since Rv ⊆ Rv(v, e), Re ⊆ Re(v, e), and Rv(v, e) ∩ Re(v, e) = ∅ for every node v and every non-incident edge e, no node can cross an edge when moving.

The property also guarantees that, given two edges, no crossings can be incurred or undone, as both cases require a node to cross an edge. Finally, Rv(v, e) and Re(v, e) can always be computed, so long as v and e are not initially collinear, cannot become collinear with l, by definition, when moving, and we are working in a continuous space.

Force System Cooling

In PrEd, nodes repel each other with inversely quadratic forces and edges attract incident nodes with linear forces. By design, distant, adjacent nodes have large attractive forces to reduce edge length. In our opinion, this behaviour can be contradictory to the aims of PrEd, as the extremities of a long edge may be tightly constrained, resulting in large, irreducible forces that delay convergence.

We overcome this problem by modifying the exponents of the force system in a way similar to previous work in revealing graph cluster structure [START_REF] Van Ham | Interactive visualization of small world graphs[END_REF]. Consider an optimal distance of δ between graph elements. Increasing the exponent of the repulsive force strengthens the repulsion of elements closer than δ and decreases the effect of elements farther than δ. Decreasing the exponent of the attractive force decreases the effect of distance, reducing the attraction of distant elements. These modifications de-emphasise the influence of distant elements.

Empirically, the best results were obtained with a gradual transition from global to local that was dependent on the iteration. Initially, we would like global forces to drive the diagram roughly towards its final shape. On the other hand, during later stages of the computation, local forces should be preferred in order to optimise the final position of elements. In the current implementation, the exponents of F r and F e (see Section 2.5) range from 2 to 4, and the exponent of F a ranges from 1 to 0.4. All exponents are increased/decreased linearly at each iteration.

As increasing exponents can generate forces of a high magnitude, we define a global distance d that cannot be exceeded when moving. This distance ranges from [3δ, 0) and is linearly decreased at each iteration. The effect of these changes increase the stability of the force system and the reliability of the parameters. They also contribute to a more even distribution of graph elements (see Figure 6).

ImPrEd New Features

In ImPrEd, several new features extend the applicability of the algorithm. With these modifications, ImPrEd can deal with a wider range of input with greater output control.

Extended Input Class

PrEd was designed to work on undirected, simple graphs. We also believe that PrEd works on connected graphs only, as some configurations of disconnected graphs can produce results that use space inefficiently. In fact, the force system of PrEd might push disconnected components apart indefinitely as there is only a repulsive force between these elements. Inspired by other force-directed algorithms [START_REF] Frick | A fast adaptive layout algorithm for undirected graphs[END_REF], ImPrEd has a gravity force F g that attracts nodes to the barycentre of the graph with constant magnitude. This force prevents the previously described behaviour and helps maintain good aspect [START_REF] Nishizeki | Planar Graph Drawing[END_REF].

ImPrEd can also multigraphs with polyline

In these the input graph converted into a graph by edge The algorithm then applied, and the result is back to its input format.

Crossable and Uncrossable Edges

In ImPrEd, each edge of the input graph can be marked as uncrossable or crossable, influencing the computation of Mv in Algorithm 1. These labels define the edge as a PrEd edge or as a crossable edge. Crossable edges contribute to the attractive force F a only. We found these edges useful for embedding graphs inside Euler diagrams. Crossable edges are temporarily removed when computing S, excluding them from M and F e computations.

Flexible and Rigid Edges

The edges λ can also be labeled as rigid or flexible. These labels define if the edges can increase or decrease their number of bends. This property is useful when edges enclose elements of the graph: by marking these edges as flexible, they can expand or contract according to the stress applied to them by their containing elements.

During the iterations of Algorithm 1, each flexible edge is polled at regular intervals to verify whether it should be expanded or contracted. Let λ be a flexible edge mapped into the path v 1 , e 1 , v 2 , e 2 . . . vn of nodes v i ∈ V and edges e i ∈ E.

When contracting λ, we look for two nodes on the path, v i , v i+2 for some i, closer than a distance α. If two such nodes exist and nodes are not inside the triangle v i , v i+1 , v i+2 , we remove v i+1 and connect v i to v i+2 . When expanding λ, we check if there are edges longer than a threshold β. If such an edge exists, we split the edge into two segments and insert them into the path. This operation is also safe, as it does not change the shape of λ. Examples of contraction and expansion of flexible boundaries are shown in Figure 7.

To obtain the best results with flexible edges, we remove node repulsion between consecutive nodes on a path. As a consequence, the attractive force exerted by an edge does not reach equilibrium and acts like a tension force that helps remove unnecessary bends. Empirically, we determined that α = 2 • δ and β = 3 • δ are good threshold values, where δ is the optimal distance between two nodes. However, during early iterations of the algorithm, this value can cause a sharp increase in the number introduced, the graph likely expanding under the effect of low attractive and high repulsive forces (see Section 4.4). This increase in the number of bends essentially increases the size of the data set and can affect the execution time of the algorithm substantially.

As an execution time/diagram quality trade-off, we could clamp the number of bends introduced into the diagram or adjust the value of α and β. For our experiments that use flexible edges, we do not limit the number of bends introduced, providing diagrams of maximum quality.

Also, it is important to note that edges should not be labeled as both flexible and crossable simultaneously. As nodes near a flexible edge often introduce bends, if the edge is also crossable, these bends would be incurred with little benefit. Therefore, in Section 6, none of our edges are simultaneously flexible and crossable. A similar argument can be made for flexible edges that initially cross, which should be avoided as well.

Results

In this section, we compare the drawing quality and execution time of PrEd and ImPrEd, on both synthetic and real data. ImPrEd (R) and ImPrEd (F) are applications of Im-PrEd with all edges marked as rigid or flexible respectively.

We ran PrEd and ImPrEd on randomly generated plane and non-plane graphs in order to test execution time stability. The plane graph data set consists of 10 plane graphs with 50, 75, 100, 150, and 200 nodes and respectively 144, 204, 294, 444, 594 edges. The non-plane data set consists of the previous data sets with 20 randomly added crossings. In all, we test 100 randomly generated graphs.

We also compared PrEd and ImPrEd on real world data. In Simonetto et al. 1). On the smallest graphs, it can count for only 5% to 10% of the total running time. In the larger diagrams, it counts for only 1% of the total running time.

Drawing Quality

Figure 9 shows the results of PrEd and ImPrEd (R) on Bertault's example graph [Ber00, Figure 1]. ImPrEd seems to reach a more stable configuration with the same number of iterations, possibly because of the new movement rules.

Figure 11 shows a random plane graph of 100 nodes optimised by PrEd, ImPrEd (R), and ImPrEd (F). Although the force cooling in ImPrEd (R) (Figure 11b) makes better use of space when compared to PrEd (Figure 11a), the drawings are very similar. The high connectivity of the graph restricts the improvements that can be made with straightline edges. ImPrEd (F) (Figure 11c) inserts bends, which increase the angular resolution of high degree nodes.

Figure 12 shows the results of PrEd, ImPrEd (R), and ImPrEd (F) on iGraphA and gGraphA. PrEd (Figures 12a and12d) produces drawings where the nodes seem forced far from the centre. In particular, Figure 12d illustrates that nodes are distributed along region boundaries. As few attractive and many repulsive forces act on disconnected nodes, we realise this effect. The original forces of PrEd are hardly usable for graphs with a large number of disconnected nodes. The results in Simonetto et al. [START_REF] Simonetto | Fully automatic visualisation of overlapping sets[END_REF] are obtained by careful selection of algorithm parameters and by increasing the repulsive force exponent (Section 4.4). However, in many of the diagrams of this article, nodes have a tendency to accumulate on set borders.

ImPrEd improves the distribution of disconnected nodes as seen in Figures 12b and12e. The modified force system of ImPrEd enables this result. Also, the δ and γ parameters of ImPrEd drive, with higher reliability, the result towards optimal node-node and node-edge distances.

Flexible edges further improve the quality of these drawings. Figures 12c and12f show improvements in the node distribution and more regular boundaries. Moreover, unnecessary bends are removed from the diagram, as shown in the close-ups, creating smoother set borders.

Finally, Figure 10 shows the results of ImPrEd (F) on gGraphB. The thick black edges, that represent the set boundaries, have smooth shapes proportional to the number of vertices contained. The grey edges, marked as crossable, helped determine the structure of the diagram and node position. However, they did not overly deform set boundaries.

Conclusions and Future Work

We presented ImPrEd, an algorithm that improves an existing graph layout while preserving the graph map. Our test cases provide some evidence that ImPrEd improves upon PrEd in terms of both running time (surrounding edges computation and QuadTrees) and drawing quality (maximal movement computation and force cooling). Through flexible edges, ImPrEd can further improve diagram quality. However, there is an execution time/quality trade-off when a large number of bends are introduced into the drawing. ImPrEd also extends PrEd to a larger input especially portions the graph lie in bounded regions. For these applications, ImPrEd can mark edges as crossable or uncrossable and as rigid or flexible, allowing for more customisable output. Finally, we showed ImPrEd can be applied to scenarios such as improving planar drawings of graphs and Euler-like diagram generation.

As future work, we would like to investigate new application fields for ImPrEd. For instance, the algorithm could be applied to combined graph-map drawings, such as GMaps [START_REF] Gansner | GMap: Visualizing graphs and clusters as maps[END_REF]. The running time improvement could also be further studied in terms of its performance with respect to edge density. However, ImPrEd still has problems scaling to large graphs. To overcome this limitation, its algorithm complexity should be further reduced as, in worst case, it coincides with that of PrEd.

Figure 1 :

 1 Figure 1: Maximal movement in PrEd. (a) The maximal movement sectors of v. (b) The corresponding Mv, containing the radius of each sector. (c) Movement limitation during the node displacement phase. Fv belongs to A 5v , so its module will be clamped to M 5 v .

Figure 2 :Figure 3 :

 23 Figure 2: Surrounding edge computation. (a) An example of plane graph. (b) The graph is divided into connected components. These components are organised into a hierarchy via containment. (c) The surrounding edges Sv are computed as the edges of the faces in Φv = {F, H}, that contain v. The edges e / ∈ Sv are dashed. (d) The disconnected boundary of the face F is computed from the connected boundaries F, according to Figure 2b.

PrEdFigure 4 :Figure 5 :

 45 Figure 4: Maximal movement computation when ve ∈ e. The dashed sectors are unlimited. (a) The result with PrEd's rules. (b) The result with ImPrEd's rules.

Figure 6 :

 6 Figure 6: Stability to the input parameters of PrEd and Im-PrEd. (a) PrEd with δ = 5 and γ = 2. (b) PrEd with δ = 2 and γ = 5. (c) ImPrEd with δ = 5 and γ = 2. (d) ImPrEd with δ = 2 and γ = 5.

Figure 7 :

 7 Figure 7: Contraction and expansion of flexible edges. The pictures show a flexible loop containing nodes. a) An overdimensioned flexible edge. c) An under-dimensioned flexible edge. b) and d) The result of ImPrEd.

Figure 8 :

 8 Figure 8: Growth of the average computation time for single iteration, with respect to the size of the random graphs presented in Section 6. a) Random plane graphs. b) Random non-plane graphs. ImPrEd (F) has not be applied to the non-plane graphs as flexible edges should not cross.

PFigure 9 :

 9 Figure 9: Layout improvement of PrEd's example graph. Both algorithms executed 100 iterations. a) Using PrEd [Ber00, Figure 1a]. b) Using ImPrEd (R).

Figure 10 :

 10 Figure 10: Result of the application of ImPrEd on gGraphB. In this graph, black edges represent the set boundaries and have been labeled as flexible. Grey edges represent element relationships and have been labeled as crossable.

Figure 11 :

 11 Figure 11: Layout improvement of a randomly generated plane graph. All algorithms used the parameters δ = 10, γ = 10, 250 iterations. a) Using PrEd. b) Using ImPrEd (R). c) Using ImPrEd (F).

Figure 12 :

 12 Figure 12: Layout improvement of graphs involved in the generation of Euler diagrams. In the first row, iGraphA. Parameters: δ = 27, γ = 21, 250 iterations. The node diameter is set to 10 to better show the graph proportions. In the second row, gGraphA. Parameters: δ = 5, γ = 4, 250 iterations. a,d) Using PrEd. b,e) Using ImPrEd (R). c,f) Using ImPrEd (F).

Table 1 :

 1 Average running times of PrEd and ImPrEd for 100, 250 and 500 iterations, over the Euler diagrams generation graphs described in Section 6. ImPrEd (PP) is the pre-processing step (Section 4.1). We only applied ImPrEd to gGraphB as the input includes edges marked as crossable, which cannot be handled by PrEd. Green numbers indicate improvements, red numbers indicate slower execution.of the final diagram. The resulting graph is used to build gGraph, which is a first approximation of the set boundaries. Finally, PrEd improves the layout of gGraph to define smoother and clearer set boundaries.Here, we consider two diagrams. The first diagram consists of 20 sets and over 2000 elements and is Figure8in Simonetto et al.[START_REF] Simonetto | Fully automatic visualisation of overlapping sets[END_REF]. Sets are films extracted from the IMDb top-chart and contain the full credited actor list.

	Graph		iGraphA	gGraphA	iGraphB gGraphB
	Nodes		149	2601	23	220
	Edges		169	507	22	362
	PrEd	100	17.50	2721.79	0.51	-
		250	40.91	6792.53	1.28	-
		500	80.46	13 585.36	2.53	-
	ImPrEd (R)	100	4.90	131.55	0.30	15.65
		250	11.81	308.83	0.73	33.49
		500	23.34	597.55	1.42	65.34
	ImPrEd (F)	100	6.68	97.22	0.47	9.25
		250	22.02	215.62	1.73	18.07
		500	35.29	409.18	3.79	30.81
	ImPrEd (PP)	-	0.13	5.64	0.02	0.14
	ImPrEd (R)	Gain	3.49	21.81	1.75	-
	ImPrEd (F)	Gain	2.25	30.90	0.73	-

  Table 1 reports the running times of PrEd and ImPrEd on the Euler diagram graphs. Again, ImPrEd (R) performs up to 20 times faster than PrEd. However, this time Im-PrEd (F) outperforms ImPrEd (R) and is up to 30 times faster than PrEd. Most probably, this is due to the high number of initial bends in the diagram that ImPrEd (F) is able to remove, effectively reducing data set size.

	On both the random data sets and Euler diagram
	graphs, pre-processing time is, for the most part, negligible
	(Table

6.1. Execution Time

Figure

8

plots running time vs. number of nodes for random planar and non-planar graphs respectively. ImPrEd (R) is consistently faster than PrEd on both data sets. On the smallest graphs, 50 nodes, ImPrEd (R) is 5 times faster than PrEd. On the largest graphs, 200 nodes, ImPrEd (R) more 20 times On other hand, ImPrEd (F) is consistently twice as slow as PrEd for all graph sizes with only small variations. As these random graphs, on average, have a high node to edge ratio, ImPrEd (F) must insert a high number of bends, slowing it considerably.

submitted to Eurographics / IEEE Symposium on Visualization 2011 (EuroVis 2011)

Acknowledgments

The second author would like to acknowledge support from Science Foundation Ireland Grant No. 08/SRC/I140.