A Schumm

M Yao

P Duvauchelle

V Kaftandjian

A Peterzol-Parmentier

X-ray imaging plate performance investigation based on a Monte Carlo simulation tool

Keywords: Monte Carlo simulation, Computed Radiography, Imaging Plate Performance, High energy applications

Introduction

Computed radiography (CR) system, based on the use of an imaging plate (IP) composed of photostimulable phosphor (PSP) powder embedded in organic binder, is one of the efficient digital replacement techniques of the traditional conventional film radiography [1,[START_REF] Rowlands | The physics of computed radiography[END_REF]. Aside from keeping the advantage of the traditional radiography such as flexibility of the detector, CR brings other advantages, which are direct digital readout image, reusability of the detector, elimination of the developing process with chemical product, reduction in exposure (through great sensitivity of IP) and dynamic range larger than 5 orders of magnitude in X-ray dose [START_REF] Seggern | Photostimulable X-ray storage phosphors: a review of present understanding[END_REF][START_REF] Miyahara | A new type of X-ray area detector utilizing laser stimulated luminescence[END_REF][START_REF] Thoms | The quantum efficiency of radiographic imaging with image plates[END_REF][START_REF] Thoms | The dynamic range of X-ray imaging with image plates[END_REF]. CR proves good performance within low energy range, but regarding the high energy range (up to MeV), which is also normal working energy range for industrial inspection of high attenuation objects, CR image quality is poor. The reason is that IP response is poor to the direct beam which carries the sharp geometry information of the object, but sensitive to secondary beam emerging from the object which is responsible for blurring. This issue has already been addressed for film radiography by the use of metallic screens which increase the dose received in the film thanks to secondary radiations occurring in the screen. Metallic screens should be placed in contact with the film so that emitted electrons can reach the argentic emulsion of the film. As far as IP is concerned, the use of metallic screens is also requested at high energies by current standards, so that the high energy photon beam could be rendered into electrons and lower energy photons, to which IP is more sensitive [START_REF]Non-destructive testing -Industrial computed radiography with storage phosphor imaging plates -Part 2: General principles for testing of metallic materials using X-rays and gamma rays[END_REF][START_REF]Non-destructive testing of welds -Radiographic testing -Part 2: X-and gamma-ray techniques with digital detectors[END_REF]. However, the commercial IPs are usually comprised of a protective layer coated on phosphor layer to provide protection against normal handling, and this layer might prevent the rendered photon or electron from entering the phosphor layer. As regards the spatial resolution, the metallic screens also further spread the incident beam, which means that a compromise should be found between intensifying the signal (by increasing the screen thickness) and improving spatial resolution. Moreover, when compared with film radiography, the spatial resolution is further degraded due to the light scattering at the phosphor grains during readout process [START_REF] Rowlands | The physics of computed radiography[END_REF][START_REF] Seggern | Photostimulable X-ray storage phosphors: a review of present understanding[END_REF]. In order to follow the interactions involved in the metallic screen and IP detector system, a Monte Carlo simulation is of particular interest. It can give insight on physical phenomena which is difficult or impossible to observe experimentally, such as scattering effect, spatial distribution of deposited energy, K absorption and K fluorescence reabsorption. Monte Carlo simulation has already been applied for system performance study and optimization as well as experimental result prediction in a number of studies [START_REF] Rogers | Fifty years of Monte Carlo simulations for medical physics[END_REF][START_REF] Boone | A Monte Carlo study of X-ray fluorescence in X-ray detectors[END_REF][START_REF] Correa | Computed radiography simulation using the Monte Carlo code MCNPX[END_REF][START_REF] Liaparinos | Monte Carlo study on the imaging performance of powder Lu 2 S iO 5 : Ce phosphor screens under X-ray excitation: Comparison with Gd 2 O 2 S : T b screens[END_REF][START_REF] Kausch | Monte Carlo simulations of the imaging performance of metal plate/phosphor screens used in radiotherapy[END_REF][START_REF] Jaffray | Monte Carlo studies of X-ray energy absorption and quantum noise in megavoltage transmission radiography[END_REF][START_REF] Liaparinos | Modeling granular phosphor screens by Monte Carlo methods[END_REF][START_REF] Liaparinos | The Monte Carlo evaluation of noise and resolution properties of granular phosphor screens[END_REF][START_REF] Souza | Methodology for digital radiography simulation using the Monte Carlo code MCNPX for industrial applications, Applied radiation and isotopes: including data, instrumentation and methods for use in agriculture[END_REF][START_REF] Li | Dose response of BaFBrI : Eu 2+ storage phosphor plates exposed to megavoltage photon beams[END_REF][START_REF] Cho | Monte Carlo analyses of X-ray absorption, noise and detective quantum efficiency considering therapeutic x-ray spectrum in portal imaging detector[END_REF][START_REF] Wu | Optimal steel thickness combined with computed radiography for portal imaging of nasopharyngeal cancer patients[END_REF], but lots of these works were conducted for medical application. Though in some works [START_REF] Souza | Methodology for digital radiography simulation using the Monte Carlo code MCNPX for industrial applications, Applied radiation and isotopes: including data, instrumentation and methods for use in agriculture[END_REF], [START_REF] Li | Dose response of BaFBrI : Eu 2+ storage phosphor plates exposed to megavoltage photon beams[END_REF] and [START_REF] Wu | Optimal steel thickness combined with computed radiography for portal imaging of nasopharyngeal cancer patients[END_REF], the spectral response of IP material and the spatial response to certain irradiation sources have been reported, there is still a lack of the detailed spatial response characterization up to high energies ( MeV). Moreover, in all previous studies, the IP system is studied alone, without taking into account the object under test. At high energies, the object is responsible for a huge amount of secondary radiation, which is not a useful part of signal as regards defect detection, and thus, image quality optimization consists in not only increasing the amount of deposited energy due to the direct beam (use f ul signal), but also decreasing the deposited energy due to the secondary beam. Direct and secondary beams must here be understood at the object level. For the purpose of understanding how the relevant operating parameters affect the X-ray image and optimizing the image quality supplied by CR system, we have developed a Monte Carlo simulation tool based on PENELOPE [START_REF] Penelope ; | Organització de Cooperació i Desenvolupament Econòmic, PENELOPE 2008: A code system for Monte Carlo simulation of electron and photon transport : workshop proceedings[END_REF]. This tool tracks separately the direct and secondary as well as photon and electron radiations. The present paper will show the use of this tool through an optimization study of a realistic inspection case. The optimization was carried out by first characterizing the response of different combinations of IP and screens. The detector response has been analyzed in terms of performance metrics such as energy absorption efficiency (EAE) and modulation transfer function (MTF) for various incident energies up to 1.3 MeV.

Simulation tool and performance metrics

Simulation method

Our simulation tool is based on the PENELOPE subroutine package, and involves a classical situation of radiography with a radiation source, an object and a detector. Here, by detector we mean the combination of IP and metallic screens since the purpose of these screens is to optimize the IP response. All physical effects are taken into account in the simulation: photoelectric effect, Compton and Rayleigh scattering, pair production, as far as photons are concerns. For electrons and positrons, Bremsstrahlung emission, elastic and inelastic scattering as well as positron annihilation are considered. The geometry adopted is shown in Figure 1(a), as well as the terms employed for the different beams issued from the object under test. After the object, the solid arrow represents the direct radiation, i.e., the part of photons which did not undergo any interaction, and all radiation produced from this beam in the following is called the use f ul beam. On the contrary, the dashed arrow represents the secondary radiation issued through the object, and all further radiations produced are denoted as not use f ul. The information is subdivided into the signals due to the photons in the use f ul beam or not use f ul beam (denoted as UPS and NUPS), the electrons in the use f ul beam or not use f ul beam (denoted as UES and NUES), and, if any, the positrons in the use f ul beam or not use f ul beam (denoted as UPoS and NUPoS). All the different particles (photons/electrons/positrons respectively use f ul/not use f ul) are tracked and registered separately during the Monte Carlo simulation process. For each of these particles, the simulation monitors two kinds of signals: the deposited energy distribution inside the IP and the step-by-step transmitted spectra after each element between source and IP (Figure 1(b)). A uniform three-dimensional orthogonal grid is used to sample the IP in such a way that the energy deposition distribution is recorded in 3D. For this reason, the IP geometry should be cuboid. The incident beam could be either diverging or not, emitted from a point or an extended source, and its spectrum can be specified (such as that emitted from a cobalt or iridium source) or be a fixed energy varying within a specified range {E min ;E max } with an increment of E incmt . The object and metallic screens could be of arbitrary material and geometry. In all the following, the IP was modeled as a multi-layered structure which is comprised of, in sequence, a 6 µm protective layer, a 150 µm phosphor layer, a 254 µm support layer and a 25.4 µm backing layer. The materials of these four layers are respectively Mylar for the protective and support layers, BaFBr : Eu 2+ with a packing factor of 60% for the phosphor layer, and polycarbonate for the backing layer. As the phosphor layer is the actual effective layer for information storing, the orthogonal grid was applied to the phosphor layer for energy deposition distribution recording.

Performance metrics

The following outputs are available from our simulation tool:

• 3D deposited energy due to each kind of beam and each kind of particle:

-Useful photon signal DEP UPS (x, y, z);

-Not useful photon signal DEP NUPS (x, y, z);

-Useful electron signal DEP UES (x, y, z);

-Not useful electron signal DEP NUES (x, y, z);

• Transmitted photon spectra after different elements (i.e. object, metallic screens and IP protective layer);

• Transmitted electron spectra after different elements.

From these outputs, the following metrics are used to investigate the IP response under different simulation configurations:

• Energy absorption efficiency (EAE). It is the fraction of X-ray energy absorbed in IP, and is calculated as the total energy deposited in IP divided by the incident beam energy. Here, the total energy is considered, i.e. the sum of all signals due to all particles. If needed, use f ul versus not use f ul energy can be separated (see the results section).

• Fraction of deposited energy due to electron beam ( f ElDep ). The f ElDep is calculated as the ratio of the deposited energy due to electrons into IP over the total deposited energy in IP. Note that the electron beam is produced by the metallic screens.

• Modulation Transfer Function (MTF) and spatial frequency at 20% MTF ( f MT F20 ). The MTF is a basic metric to express the system spatial resolution. It is calculated as the modulus of the Fourier transform of the Point Spread Function (PSF) normalized to unity at zero spatial frequency. In practice, the Line Spread Function, describing the system response to a line, is usually measured, and then the corresponding one-dimensional MTF is [START_REF] Liaparinos | Monte Carlo study on the imaging performance of powder Lu 2 S iO 5 : Ce phosphor screens under X-ray excitation: Comparison with Gd 2 O 2 S : T b screens[END_REF]:

MT F( f ) = F (LS F(x)) = F         y PS F(x, y)         = F         z,y DEP(x, y, z)         (1)
where f is the spatial frequency, and F denotes Fourier transform, and DEP(x, y, z) here is obtained with a point The geometry adopted for our simulation is a classical radiography chain (not drawn to scale). A Cartesian coordinate system is used to track the particle transport. The y-axis points outward. A uniform three-dimensional orthogonal grid is applied to IP to record the deposited energy distribution. (b) Simulation outputs: 3D map of deposited energy in IP DEP(x, y, z) and an example of spectrum after one of the elements between source and IP.

non-diverging source (pencil beam). The f MT F20 is the spatial frequency at which the MTF equals 20%. Depending on the investigation required, it is possible to consider the energy due to the different combinations of use f ul/not use f ul and photon/electron.

• Energy absorption distribution along IP depth (EAD z ) and its centroid (EADC). An originality of our tool is to consider the distribution of deposited energy along the depth of IP (z direction). This is not common in other tools. Our idea is following: the latent image stored in IP (i.e. the deposited energy) is read by a laser light, and while the laser light enters IP, it diffuses along its traveling path due to the absorption, reflection and scattering effects, which contributes to image blurring. The laser power used during readout modifies the penetration of laser light, and thus, has an influence on final spatial resolution. Depending on the deposited energy profile along IP depth, an optimal laser power could be determined. For this reason, we would like to study the energy deposition distribution along the IP depth.

The EAD z is a function of z and is computed as the sum of DEP(x, y, z) over x and y axes normalized by the incident beam energy. The centroid of EAD z could be understood as the average position of the deposited energy, which is calculated using the following formula:

EADC = z z • EAD z (z) z EAD z (z) (2)

Context of our study: typical inspection case

Our typical industrial inspection case consists of irradiating a thick object (70 mm steel plate) with a high energy γ source (Cobalt). In order to explain the image quality issue, we show the result obtained when IP is used alone as a detector, without any metallic screen to intensify the signal. A pencil beam is used to irradiate the object, in order to show the spread out responsible for spatial resolution degradation. IP is placed directly behind the object, as in the industrial case. Figure 2 illustrates the deposited energy images in XZ plane due to respectively the use f ul photon signal (UPS), use f ul electron signal (UES), not use f ul photon signal (NUPS) and not use f ul electron signal (NUES). Remember that z direction is the incident beam direction, and is limited to 150 µm, which is the thickness of the phosphor layer. As can be seen, the use f ul signal images from photons and electrons are quite sharp and clear, while the not use f ul signal images are blurry. In order to compare the amount of energy, in Figure 3, we plot the deposited energy along z direction for the respective signals.

It is worth noting that the not use f ul signals are dominant, as only about 15.7% of deposited energy is due to use f ul signals. When looking at the depth penetration, we see that the electrons deposit their energy preferably towards the front side of IP (z = 0), while the photons near the rear side (z=150 µm). In Figure 4, we investigate the contribution of the different signals to the overall MTF curve (dash-dotted line with point marker). Although the use f ul signals UPS and UES have good MTFs (solid lines), the overall MTF is predominantly influ- enced by the not use f ul signals (NUPS, NUES). This is the quantitative aspect of the visual spread out already shown in Figure 2. It must be noted here that in real life, electrons might not have as much influence as what is shown here, because air was not considered in the simulation. However, this does not change the conclusion that not use f ul photons degrade a lot spatial resolution. Thus, to conclude on our industrial context, we need a detector which should be efficient to the high energy use f ul signal, but insensitive to the lower energy not use f ul signal. The exact role of electrons must be very carefully tracked, as they might carry quite a lot of energy, but could be absorbed in the protective layer of IP itself. Thus, we see that a fine detector response characterization must be performed. A detailed analysis on the spectral and spatial response will be presented in the next section for physical phenomena comprehension and determination of optimal detector configuration.

Detector response characterization

We now investigate the impulse response of the system, which means the response to a pencil beam. All the elements (metallic screens and IP) were set cuboid, and no object is included in order to characterize only the detector system. The setting of all parameters is listed in Table 1. The incident number of photons is 10 7 for all the simulations considered in this section. Two combinations of screens were compared with IP alone. The first screen configuration (lead screen) is that requested by current standards. The second configuration (lead and tin) is chosen as an attempt to absorb K fluorescence of lead.

Spectral response analysis 4.1.1. Investigation of transmitted step-by-step spectra

In order to investigate the physical phenomena involved in each screen, we follow the photon spectrum after each step of the complete combination of screens (lead and tin). A monochromatic incident beam of 150 keV is chosen in such a way to be in the energy range where fluorescence occurs. Figure 5 illustrates the photon spectra:

• The first step in the detector system is the lead layer (asterisk) where it is visible that the initial monochromatic beam results in a continuous spectrum of secondary scattered radiation and fluorescence lines (K lines at 87, 84 and 74, 72 keV and L lines around 15 keV).

• The second step is the tin screen (dashed line): first of all it is observed that tin can effectively absorb the lead fluorescence (about 2 orders of magnitude) and also a big part of the scattered radiation between 30 and 80 keV. However, tin also yields its own fluorescence lines at around 30 keV.

• The last step corresponds to the protective layer of IP (plus sign): the overlapping of the dashed and cross marker curves indicates that no effect is introduced by the protective layer on the photon spectrum. Figure 6 shows the electron spectra after the different layers. In the photoelectric effect, an excitation photon of energy E is absorbed by the atom, which leads an excited atom and a photoelectron leaving the atom with a kinetic energy E k = E -U i , where U i is the ionization energy of the electron [START_REF]Hand book of medical imaging[END_REF].

• The first step (asterisk) represents the electron spectrum after the lead screen, where two peaks are visible, respectively at 62 keV and 135 keV, corresponding to the photoelectrons coming from the lead K-shell (whose ionization energy is about 88 keV) and L-shell (whose ionization energy is about 15 keV).

• After the tin layer (whose K-and L-shell ionization energies are about 29 keV and 4 keV respectively), the photoelectron peaks are observed at 121 keV and 146 keV (dashed line).

• When comparing the spectra after tin (dashed line) and after the protective layer (plus sign), we can observe that the photoelectron peaks are much reduced, but overall, the energy stopped by the protective layer is negligible with respect to the total incident energy into this layer. Thus, in this energy range, the protective layer do not prevent electrons from entering the phosphor layer.

The study of the spectra is illustrative and helps to understand phenomena, but for high energies, a more global study of absorption efficiency is needed.

Energy absorption efficiency

Figure 7(a) shows the energy absorption efficiency plots of the 3 different detector configurations.

• The response of IP alone (solid line) is quite efficient, nearly 1, at low energies, but with the increase of excitation energy, the curve drops quickly towards zero. The small rise occurring between 20 and 30 keV corresponds to the barium absorption edge.

• The lead screen (dashed line) filters the low energy radiation but intensifies the high energies: the dashed curve is above the curve of IP alone for energies greater than 131 keV. The barium absorption edge has disappeared due to the filtering of low energies. At around 89 keV a quick drop arises which is due to the strong K-absorption of lead. • With the combination of lead and tin screens (solid line with diamond markers), similar effects are observed, except that the signal intensification starts from 226 keV, and an additional drop appears at 30 keV which corresponds to the K-absorption of tin.

In order to compare the intensification effect of the different combinations of metallic screens, the absorption curves of the two configurations were normalized by that of IP alone (Figure 7(b)). Starting from 131 and 226 keV respectively, the two kinds of screen combinations intensify the absorbed signal in IP (normalized absorption to 1), and the intensification increases continuously. Regarding only the absorption intensification in the high energy range, the detector configuration with only lead is preferred, as it has the highest absorption level among the three. However, the tin plus lead combination has a better filtering of low energies. Thus, when considering the compromise introduced in section 3, namely the need to increase high energy weight (in order to absorb the use f ul signal) and avoiding the low energy not use f ul signal, the two screen combinations can be discussed (see section 5). On the other hand, the spatial resolution must also be compared, which is the aim of the next section.

Spatial response analysis

4.2.1. MTF curve and the limit of resolution at 20% MTF ( f MT F20 ) Figure 8(a) shows the MTF curves of the IP alone for different excitation energies. In the low to medium energy range (see dark blue curves, up to about 400 keV), the MTF curve drops quickly with the increase of excitation energy, whereas, in the high energy range (see light blue to dark red curves, 400 keV to 1300 keV), the MTF curves overlap each other, meaning that the spatial resolution could be considered energy independent in this high energy range. Two gaps appear at 14 keV and 38 keV resulting from the K fluorescence reabsorption of brome and barium, as the fluorescence emission is isotropic, and thus, the PSF spreads out. Figure 8(b) is the plot of f MT F20 versus excitation energy for IP alone. The two slumps correspond again to the K fluorescence reabsorption of brome and barium (like the two gaps just described above). As the excitation energy continues to increase, exceeding the K-edge energy, less fluorescence is produced, so less K fluorescence reabsorption, and hence the limit of resolution improves right after the K-edges, which was difficult to see in Figure 8(a) due to the color scale. Figure 9 compares the resolution of the three detector configurations. With the lead screens on both sides of IP, a degradation of spatial resolution is introduced, especially at 89 keV, at which the fluorescence of the lead starts to yield. As a photon with a higher energy has a larger free traveling path, the lead fluorescence could spread further than that of barium and brome, and hence the resolution drop is much greater at 89 keV than it is at 14 and 38 keV. At energies lower than 89 keV, the addition of a tin screen between lead and IP further degrades the resolution since it is another layer spreading the incident beam. However, when the excitation energy exceeds 89 keV, as the lead fluorescence could be effectively attenuated by the tin screen (Figure 5), the resolution improves. This improvement, however, decreases as the excitation energy continues increasing, and becomes negligible at 1300 keV. As a conclusion of this part, even considering only spatial resolution, it is difficult also to choose the best screen configuration because the behavior respectively at low to medium and high energies is not the same, and both should be considered for our typical industrial case.

The energy absorption distribution along z (EAD z ) and

its centroid (EADC) Now we would like to investigate the energy deposition with respect to the IP depth (z direction).

Figure 10(a) is a 2D visual image of energy absorption (in a color scale) as a function of excitation energy in abscissa, while the IP depth is represented in the vertical direction, with the front side at z = 0 and the rear side at z = 150 µm. First of all, it is very clear that absorption is better in the low energy range (red color is only near 0 keV), and most of the energy is deposited towards the front side of IP. In order to have a better visualization of the influence of z direction, Figure 10(b) shows the profile of energy deposition along z, EAD z , at 20, 40, 100 and 1000 keV. We see that when excitation energy is small, the deposited energy is mostly near the front surface of IP, while with the increase of the excitation energy, the distribution is more and more towards the rear surface. The dropping at the boundaries is probably due to particle escaping.

In order to facilitate the comparison, we plot the centroid EADC versus excitation energy in Figure 11. In our case, the centroid value is between 0 and 150 µm, where 0 µm means that all energy is deposited at the front surface of IP, and 150 µm the rear surface. When viewing the response of IP alone, at low energies, the energy is preferentially deposited on the front side, and with the increase of excitation energy, the energy tends to be deposited towards the rear side so that the centroid value in-

creases.

When comparing screens with IP alone, we see that the centroid is shifted towards the front side to a different extent. Similarly to the absorption efficiency and sharpness measurements, it is also observed a sudden change due to the K-absorptions of fluorine, brome and barium in the solid line, and lead in the dashed line. When comparing both screens, a small decrease of the centroid is observed with only the lead screen, which must be elucidated. It must be noted that the centroid investigation is important because the laser readout process is usually done from the front side of the IP and limited to a small depth. Thus, if we show that the centroid is lower, there is an interest for a best readout process. As an attempt to elucidate the decrease in the centroid distribution with lead, in Figure 12(a), we study separately the centroid obtained for the photon and electron beams: the photon beam EADC curve overlaps the EADC curve of IP alone, whereas the electron beam EADC is closer to the front surface, and it increases with excitation energy. Here we see that the electron beam is mostly responsible for shifting the deposited energy towards the front side, but the degree of the shift depends on the fraction of the deposited energy which is due to electron beam. Thus, in Figure 12(b) we plot the fraction of the deposited energy due to the electron beam ( f ElDep ): the f ElDep also increases with excitation energy and at energies over 400 keV, more than half of the deposited energy comes from the electrons emerging from the lead screen. Thus, when comparing lead screen with IP alone, at low energies, even though the electron beam centroid is very close to zero, as the fraction of deposited energy is small, no obvious shift is observed; with the increase of excitation energy, the centroid shift could be observed, and the optimal centroid is at around 300 keV, an energy for which the electron beam centroid is still small and the fraction of deposited energy is important.

Results: industrial case optimization

According to the detailed analysis in the previous part, we find that both the two detector configurations with metallic screens could increase the use f ul signal weight in the final image since they could effectively increase the absorption efficiency at high energies and eliminate the influence of low energy signal, especially IP + 0.2Pb0.8S n which has better filtering effect at low energies. While limiting the blurring introduced by not use f ul signal, the screens also degrade the spatial resolution, therefore, a suitable detector should provide good use f ul to not use f ul signal ratio without losing much spatial resolution. As experiments, the outcome of a Monte Carlo code is affected by statistical uncertainty, which is directly related to the number of photons involved in the process. The greater the number of incident photons, the less the relative uncertainty. The error is usually estimated as three times the standard deviation of the deposited energy. In our case, the number of initial photons has been chosen in order that this error is relatively small. In Figure 13 we show the use f ul over not use f ul signal ratio: the ratio is greatly improved with IP + 0.2Pb0.8S n. The incident photon numbers of the three configurations are respectively 6 × 10 8 , 3 × 10 8 and 3 × 10 8 . In Figure 14 we see that the MTF with IP + 0.2Pb ( Figure 14(a)) and with IP + 0.2Pb0.8S n (Figure 14(b)) are nearly the same. Comparing the spatial frequency at 20% MTF, IP + 0.2Pb0.8S n (1.62 lp/mm) is a little better than IP + 0.2Pb (1.45 lp/mm). Hence IP + 0.2Pb0.8S n is overall better for our industrial case.

Discussion

Through the IP response characterization and a detailed analysis on the spectral and spatial responses, we optimized an inspection case. IP+0.2Pb0.8S n is the overall preferable detector configuration since it increases the most the high energy use f ul signal weight without losing too much the energy absorption efficiency and spatial resolution. From the performed study, we find that the metallic screens could: enhance the high energy signal and filter the low energy signal, degrade spatial resolution (the fluorescence of the screens is crucial), and shift the deposited energy toward the front surface of IP. The electron beam is mostly responsible for shifting the energy deposition distribution, and the degree of the shift depends on both the fraction and the centroid of the deposited energy due to electron produced by the metallic screens.

Conclusion

In this paper, we presented a Monte Carlo simulation tool for investigating the physical phenomena involved in the energy deposition process in IP. This tool allows different case studies on understanding the experimental observations and optimizing the image quality and choose the best IP/screen combination. This tool also provides the evaluation of the spectra transmitted along the different elements so that the fluorescence and scattered radiation production could be monitored as well as the influence of the protective layer of the IP. As a future work, a comparison between experimental and simulation results of the detector performance (i.e. spectral and spatial response) will be of great interest. The IP response characterization results obtained by this simulation tool will be applied to our future work on CR imaging chain modeling and optimization.

Figure 1 :

 1 Figure1: (a) The geometry adopted for our simulation is a classical radiography chain (not drawn to scale). A Cartesian coordinate system is used to track the particle transport. The y-axis points outward. A uniform three-dimensional orthogonal grid is applied to IP to record the deposited energy distribution. (b) Simulation outputs: 3D map of deposited energy in IP DEP(x, y, z) and an example of spectrum after one of the elements between source and IP.

Figure 2 :Figure 3 :

 23 Figure 2: Deposited energy image in XZ plane, due to UPS, UES, NUPS and NUES.

Figure 4 :

 4 Figure 4: MTFs of the deposited energy image due to UPS, UES, NUPS and NUES, and the MTF of total deposited energy image.

  2, 3: Pb L lines Peaks 4, 5, 6, 7: Pb K lines Peaks 8, 9: Sn K lines

Figure 5 :

 5 Figure 5: Step-by-step photon spectra obtained with Monte Carlo simulation.

Figure 6 :

 6 Figure 6: Step-by-step electron spectra obtained with Monte Carlo simulation.

Figure 7 :

 7 Figure 7: (a) The photon absorption efficiency as a function of energy for different detector configurations and (b) the absorption efficiency curves of different detector configurations normalized by that of IP alone.

Figure 8 :Figure 9 :

 89 Figure 8: (a) MTF curves at different excitation energies for IP alone and (b) the f MT F20 versus excitation energy. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Figure 10 :Figure 11 :Figure 12 :

 101112 Figure 10: (a) 2D visual of absorption versus IP depth and excitation energy and (b) the EAD z at 20, 40, 100 and 1000 keV. 0 µm is the front side of IP and 150 µm the rear side. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Figure 13 :Figure 14 :

 1314 Figure 13: The ratio of use f ul over not use f ul signal.

Table 1 :

 1 Parameter setting. The IP was between the front and back screens, and the screens were placed in contact with IP

			Elements	IP alone	Configurations IP+0.2Pb	IP+0.2Pb0.8Sn
	Source	Type:photon source Energy range in keV (E min , E max , E incmt )	fixed energy (2,1300,1)	fixed energy (2,1300,1)	fixed energy (2,1300,1)
	Object	Object	none	none	none
		Front screens	none	0.2mmPb	0.2mmPb +0.8mmSn
	Detector	IP	NbBX×NbBY×NbBZ 100×100×100 100×100×100 100×100×100 Voxel dimension (µm 3 ) 10×10×1.5 10×10×1.5 10×10×1.5
		Back screens	none	0.2mmPb	0.8mmSn +0.2mmPb

Preprint submitted to Spectrochimica Acta, Part B November 20, 2014


Florian Montagner

Valérie Kaftandjian

Delphine Farlay

Daniel Brau

Georges Boivin

PhD Hélène Follet
email: helene.follet@inserm.fr

Validation of a novel microradiography device for characterization of bone ineralization

Keywords: mineralization of bone, quantitative microradiography, digital detector, spatial resolution, contrast to noise ratio

de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Bone is a living and composite material. Its intrinsic changes affect its ability to resist fracture. Bone is composed of organic matrix (~30%, mainly type I collagen), mineral (~60%) and water (~10%) [START_REF] Farlay | Bone Mineral Quality[END_REF]. The inorganic composition of bone ("bone mineral") is a poorly crystallized apatite which slightly differs from stoichiometric hydroxyapatite (HA) Ca10(PO4)6(OH)2 [START_REF] Blanchard | Intravoxel bone micromechanics for microCT-based finite element simulations[END_REF]. The density of the mineral part of bone is an important determinant of bone strength [START_REF] Follet | The degree of mineralization is a determinant of bone strength: a study on human calcanei[END_REF][START_REF] Boivin | Alendronate increases bone strength by increasing the mean degree of mineralization of bone tissue in osteoporotic women[END_REF]. Bone mineral quality can be measured at the organ (i.e., bone and surrounding tissues) or tissue levels (i.e., mineral part of bone only) with different methods. At the organ level, the current and clinical reference method for measuring the amount of mineral in bone is by dual energy X-ray absorptiometry (DXA) which gives the so-called bone mineral density area (aBMD) [START_REF] Monija | Combined Vertebral Assessment and Bone Densitometry Increases the Prevalence and Severity of Osteoporosis in Patients Referred to DXA Scanning[END_REF][START_REF] Blake | DXA in Adults and Children[END_REF]. The term "area" comes from the fact that DXA is a two-dimensional method and thus only gives the bone mineral density per unit area (g.cm -2 ). A three-dimensional measurement of the mineral density can be obtained by peripheral quantitative computed tomography (pQCT) [START_REF] Glüer | Quantitative Computed Tomography in Children and Adults[END_REF][START_REF] Engelke | Quantitative computed tomography (QCT) of the forearm using general purpose spiral whole-body CT scanners: Accuracy, precision and comparison with dual-energy X-ray absorptiometry (DXA)[END_REF]. As a volumetric method, this measurement is called vBMD (for "bone mineral density volume") and is given in g.cm -3 . This measurement concerns only the zones of radius and tibia. At the tissue level, the main issue is to avoid any geometrical effect. Measuring only the mineral part is defined by the tissue mineral density parameter (TMD in g.cm -3 ), which avoids the porosity, the vascular and canalicular network. TMD is obtained by high resolution pQCT and microcomputed tomography (µCT) but considers only the bone mineral part. In addition to other methods for measuring the mineral part of bone [START_REF] Boivin | Methodological considerations in measurement of bone mineral content[END_REF], quantitative microradiography allows to obtain the so-called degree of mineralization of bone (DMB in g.cm -3 ), defined as the density of HA [START_REF] Jowsey | Quantitative Microradiographic Studies of Normal and Osteoporotic Bone[END_REF]. This well known and referenced technique is based on the radiographic acquisition of an embedded bone section of known thickness [START_REF] Boivin | The degree of mineralization of bone tissue measured by computerized quantitative contact microradiography[END_REF]. Since the method uses high resolution radiographic films, the spatial resolution is very high and the DMB obtained is very accurate. However, due to the disappearance of manufactured radiographic films, this technique needs to be updated using digital sensors. The purpose of this study was to setup the upgrade of microradiography from high resolution radiographic films to a digital detector, and to quantify the image quality (in terms of contrast, noise and spatial resolution). After presenting the principle of DMB measurement, the specifications of the two systems are shown, as well as the corresponding calibration curves. Then the image quality obtained with the two systems was investigated to validate the replacement of the film by a digital detector.

Theoretical background for DMB measurement

When a monochromatic X-ray beam of N0 photons of energy E irradiates a sample, the number of transmitted photons N(E) follows the Beer-Lambert law (Eq. ( 1)).

    d ρ ρ Z E µ e E N E N      ) , , ( 0  (1)
where µ/ρ is the mass attenuation coefficient of the sample (cm 2 .g -1 ), Z its atomic number, ρ its density (g.cm -3 ), and d its thickness (cm). The Beer-Lambert law indicates that the global attenuation of an X-ray beam depends on its energy and material-dependent parameters, among which some are independent of energy (thickness and density), and one is energy dependent (µ/ρ). The mass attenuation coefficient µ/ρ is introduced in the equation instead of the linear attenuation coefficient µ because it is independent from the density of the sample and, thus, on its physical state. µ/ρ values are well known for all simple chemical elements and can be computed for any compound, provided that its chemical composition is known.

The degree of mineralization of bone (DMB) is defined by the density of HA ρHA [START_REF] Rowland | Microscopic metabolism of calcium in bone. III. Microradiographic measurements of mineral density[END_REF] and computed from the number of photons transmitted by a bone sample. With the incident flux being unknown, the measurement is based on a comparison with a reference object. Aluminum (Al) was chosen as a reference because it is easy to use and perfectly known (thickness and density) [START_REF] Wallgren | Biophysical analyses of the formation and structure of human fetal bone; a microradiographic and x-ray crystallographic study[END_REF].

The Beer-Lambert law applied respectively to Al and HA gives the following Eqs (2) and (3):

    AL AL AL d ρ ρ Z E µ AL e E N E N              ) , , ( 0  (2)     HA HA HA d Z E µ HA e E N E N                 ) , , ( 0 (3)
Eqs ( 2) and [START_REF] Follet | The degree of mineralization is a determinant of bone strength: a study on human calcanei[END_REF] show that the number of photons transmitted by the Al or bone sample is related to thickness and density. In our case, the unknown parameter is bone mineral density (ρHA), while all others parameters are known. Thus, if the signals measured on bone and Al are identical, which means that NAL=NHA, considering the same incident flux N0, we get Eq. ( 4):

AL AL AL HA HA HA d Z E µ d Z E µ                            ) , , ( ) , , ( (4)
This allows us to link the DMB (ρHA) to the Al thickness dAL (5):

AL HA AL AL HA AL HA d K d d Z E µ Z E µ DMB . ) , , ( ) , , (   
                        (5)
To apply Eq. ( 5), the K factor must be known. The density of Al is perfectly known (2.7 g/cm 3 ), as well as the sample thickness dHA (100 µm). The mass attenuation coefficients of Al and HA depend on energy, and their values are tabulated [START_REF] Hubbell | Tables of X-ray mass attenuation coefficients and mass energy-absorption coefficients 1 keV to 20 MeV for elements Z=1 to 92 and 48 additional substances of dosimetric interest[END_REF]. The two curves of mass attenuation coefficients of Al and HA versus energy are shown in Fig. 1. For an energy of 8.05 keV, the ratio of mass attenuation coefficients is 0.571 and the K factor value is 0.571*2.7/0.01=1.54.

Specifications of the two systems using radiographic film and digital detector

As explained in the previous section, the DMB measurement depends on the energy. A monochromatic beam is preferred in order to compute the K factor (Eq. ( 5)). In practice, for the two systems, an X-ray tube is used with a copper anode (Kα radiation energy is 8.05 keV) and a nickel filter, like the very first studies of quantitative microradiography [START_REF] Wallgren | Biophysical analyses of the formation and structure of human fetal bone; a microradiographic and x-ray crystallographic study[END_REF][START_REF] Holmstrand | Biophysical investigations of bone transplants and bone implants; an experimental study[END_REF]. The K-edge absorption of nickel is 8.3 keV, which is just above the Kα characteristic line of copper. Thus, the nickel filter is well adapted to attenuate energies just above the copper characteristic line. The spectrum obtained after filtration is still polychromatic, but with a strong preponderance of the copper characteristic line at 8.05 keV. The two microradiographic systems are illustrated in Fig. 2. With the system using radiographic films, samples are directly in contact with the film. Thus, no magnification is used and no geometric unsharpness occurs. As a consequence, in contact microradiography, the focal spot size of the X-ray tube can be large, thus allowing a high radiation intensity. On the contrary, with the digital detector, since the detector pixel size is higher than the desired spatial resolution, the magnification factor must be increased in order to obtain a smaller pixel size at the sample level. However, the geometrical unsharpness increases as well, and the use of a microfocus X-ray tube is required. The power is therefore reduced and the radiation intensity is low. A compromise needs to be found regarding image quality because a high intensity improves the signal to noise ratio (SNR), while increasing spatial resolution requires a decrease in intensity. This compromise will be discussed in section 5. The technical specifications of the two systems are summarized in Table 1.

Experimental validation using an Al step wedge

An Al step wedge was exposed using both techniques (8 regular steps from 12.5 µm to 100 µm). Typical grey-level profiles obtained by the two systems are shown in Fig. 3. Both profiles were normalized by their maximal value in order to be shown on the same scale. However, it is worth noting that the dynamic range of the detector is 4096 (12 bits), while 255 levels are available for the film system (8 bits).

For the digital detector system, the grey-level is proportional to the energy deposited in the detector, which is related to the energy transmitted by the object, weighted by the detector efficiency. Thus, a highest signal is obtained for small thicknesses where many photons are transmitted. For the other system (radiographic film coupled to a microscope), the deposited energy in the film increases its optical density (film darkening), and the obtained signal increases with Al thickness (which means when low energy is transmitted).

Fig. 3. Normalized grey levels obtained for an aluminum step wedge with both systems. M8 represents the thickest step and M1, the thinnest.

From the experimental evolution of the measured signals, we could deduce a logarithmic relation between Al thickness and grey levels, as given by Eq. ( 6):

  b GL a d Al    ln (6)
where a and b are constant coefficients. In the case of the digital detector, the coefficient a is negative (when Al thickness increases, signal decreases and grey levels become black), and it is positive in the case of radiographic film (the light intensity is measured through the radiographic film). Finally, Eqs. ( 5) and ( 6) allowed us to obtain a relation between grey levels and DMB as follows:

* ) ln( * b GL a DMB    (7)
where a*=K.a and b*=K.b. The obtained DMB as a function of grey levels is plotted in Fig. 4. This calibration curve allowed us to determine the DMBs for all possible values of grey levels in the range of interest. For the system using radiographic films, the gain of the microscope camera was chosen in order to optimize the dynamic range in the domain of grey levels which were close to the bone tissue values. This is why five steps were used for the whole grey-level range. The calibration was validated with the digital system by exposing a sample of 98% pure HA with a density of 3.08 g.cm -3 . The density calculated from the procedure was 3.015 g.cm -3 , i.e., a relative difference of 2.1%. Fig. 4. Calibration curves of DMB as a function of grey-level, with film and digital detector.

Quantification of the image quality based on the contrast to noise ratio (CNR) and spatial resolution

CNR

The DMB is calculated from the grey-level. The smallest distinguishable DMB values are directly related to the grey-level contrast resolution. The noise level must be taken into account, especially for digital systems. An essential image quality parameter is the CNR. The CNR is defined as the grey-level difference between two thicknesses of a reference sample divided by the noise level [START_REF] Duvauchelle | A computer code to simulate X-ray imaging techniques[END_REF]. CNR is given by the relation in Eq. ( 8):

2 2 2 1 2 1 GL - CNR GL GL GL     (8)
where: GLi is the average grey level in zone i σGLi the standard deviation of grey level in zone i.

A CNR of one means that the grey-level difference is equal to the noise level, i.e., the two thicknesses will be undiscernible. The higher the CNR, the better is the image quality in terms of contrast resolution. CNR increases when the contrast improves, which is known to be related to the energy of the X-ray beam: a low energy improves the contrast. On another hand, the CNR improves when the signal to noise ratio (SNR) improves, which is known to be related to the number of photons. Thus, the best quality should be obtained at a low energy and a high photon intensity. In our case, when increasing the high voltage, the mean energy of the spectrum increases, but the spectrum is still dominated by the copper characteristic line. A higher voltage increases also the total number of photons present in the spectrum, which increases the SNR.

Overall, the CNR is better when increasing the high voltage.

In our study, CNR was calculated between two Al steps of 62.5 and 75 µm, since they corresponded to the closest values of DMB generally observed. In order to optimize the settings of the X-ray tube, we performed two tests: (1) the high voltage was changed while keeping a constant power of 2 W (which means that the intensity decreased when high voltage increased, keeping the same integration time of 7 s); and (2), high voltage was changed with a constant exposure of 0.35 µA.s (which means the same intensity and integration time, but an increasing power). Fig. 5 shows that in both cases, CNR increased with high voltage, for the reason given above. As a matter of comparison, CNR obtained with radiographic film was about 6.5 with an acquisition time of 20 min at 25 kV, 25 mA (optimized setting values for film [START_REF] Boivin | The degree of mineralization of bone tissue measured by computerized quantitative contact microradiography[END_REF]).

With regard to the digital system settings, a high voltage of 40 kV was selected (just before the maximal value of 45 kV). Concerning the current, our objective was to preserve the best spatial resolution (as the focal spot size increased with current), a low power had thus to be selected (2 W, i.e., 50 µA at 40 kV).

To maximize the flux of photon, it is possible to increase the detector integration time and average several images, at the expense of total acquisition time. Indeed, the CNR increases with the square root of the number of integrated images. An average of 6 pictures was chosen here as a good compromise. All these parameters allowed to achieve a CNR of 5.5 which was acceptable, although less than the film CNR.

Spatial resolution

Spatial resolution measures the ability of a system to separate details which are close together.

It is influenced by the resolution of the detector and by the geometric unsharpness which depends on the magnification factor and the focal spot size (the latter being related to the tube intensity). Following the Nyquist sampling theorem, a digital detector is able to discern an object whose size is at least greater than twice the size of a pixel. In our case, the detector pixels have a size of 9 µm, thus the maximum intrinsic resolution (i.e. the resolution strictly due to the detector, without magnification and thus without geometric unsharpness) is expected to be 18 µm.

A mean for estimating the spatial resolution is to calculate the modulation transfer function (MTF) of the system [START_REF] Bengt | Resolution, unsharpness and MTF[END_REF]. This is obtained using a regular edge object of high attenuation, and then to realize the Fourier transformation of the derivative curve profile. MTF curves for the film and the detector with different magnification are illustrated in Fig. 6. Using the digital detector, the closest line corresponding to the film system is obtained for a magnification of 3.7.

A useful criterion to compare several systems is to assess the spatial resolution limit as the one given for a certain MTF value, typically 10% [START_REF] Kaftandjian | A comparison of the ball, wire, edge and bar/space pattern techniques for modulation transfer function measurements of linear X-ray detectors[END_REF]. For a magnification of 3.7, for instance, the spatial resolution limit is 100 lp/mm, i.e., 5 µm. At this magnification level, the object pixel size is 2.43 µm (detector pixel size divided by the magnification factor), which means a Nyquist value of 4.86, in good agreement with the result, indicating that resolution is limited by pixel size and not by the scintillator effect. The use of a calibration bar-space pattern (gold bars engraved in a plastic sheet) is another way to determine and verify the spatial resolution of the system. The contrast transfer function (CTF) is a plot of the contrast obtained versus the spatial frequency. The procedure of measurement is represented in Fig. 7: from the image of the reference object, a grey-level profile is extracted and allows computing the contrast obtained for each spatial frequency. CTF is always above the MTF because it corresponds to a rectangular pattern, while by definition the MTF corresponds to a sinusoidal one. CTF curves obtained for different magnifications with the digital system (full lines) and with radiographic film (dashed line) are shown in Fig. 7(c). This confirms that a magnification of 4 is necessary with the digital system to correspond with film.

It must be noted that the spatial frequency is limited to 20 lp/mm due to the minimal bar-space of the pattern which is 25 µm, and thus the curve cannot be plotted above this value. Given the MTF, the CTF can be also determined by the Coltman formula [START_REF] Coltman | The Specification of Imaging Properties by Response to a Sine Wave Input[END_REF] with the use of Eq. ( 9):

                           ... 11 11 9 9 7 7 5 5 3 3 4 M f M f M f M f M f M f C  (9)
where, M(f) = sine wave MTF C(f) = bar target CTF f = spatial frequency.

The CTF values are in agreement with the MTF values using this relation.

Comparison of an image obtained by the two systems on the same bone sample

Fig. 8 shows images performed with the two systems on an identical bone sample (Fig. 8(a)). For the digital system, two images were taken: one with a magnification of 3.7 (Fig. 8(c)) and another with the maximum magnification 10.8 (Fig. 8(d)). With the magnification of 3.7, the system allowed to obtain overall the same structures than with the radiographic film (Fig. 8(b)), but the image was less sharp. With the magnification of 10.8, the picture obtained was very similar to the radiographic film. For the three configurations, the grey-level profiles were plotted along two similar bone structures and they were very similar (Fig. 8(e)). CNR obtained with radiographic film and with digital detector (M=3.7 and M=10.8) were respectively 8.2, 7.1 and 7.9.

Conclusion

The digital system provides a good CNR (about 5.5 on the range referred DMB with conditions specified above). Moving the sample holder allows adjusting the magnification between 1.6 and 10.8 and this allows to achieve a similar spatial resolution compared to film for a magnification factor of 4. A higher magnification allows to achieve an even better MTF than film. Final images on bone samples show the same level of detail between film and digital detector at the highest magnification factor (10.8). On the other hand, the digital system dynamic range is 12 bits (i.e., 4096 values) versus 8 bits (255 values) for system using films. So the grey level information is more precise with the digital system. Finally, acquisition settings have been fixed as mentioned in Table 2. The use of a digital detector for radiographic imaging has the advantage of being easier to use than films. Indeed, film chemical revelation can be modified by temperature, user manipulation and fixator time settings. Moreover, there is no consumable anymore. We can also mention the problem of reciprocity of films: the efficiency of a digital detector is linear with exposition, but with the film, the efficiency decreases with exposition. This phenomenon is known as the Schwarzschild effect [START_REF] Djouguela | The Schwarzschild effect of the dosimetry film Kodak EDR 2[END_REF]. So with the digital system, we can easier control the evolution of photon flux when exposure time is modified compared to a system with radiographic films. Another point with the two systems is the fact that several images are required to cover the entire sample, especially when the magnification factor is higher. Thus, the total acquisition time will depend of number of images required. For a sample which needs 12 images, for instance, about 8 min are necessary against 20 for the entire sample with the radiographic film.

In conclusion, the two systems have similar performances with regard to image quality (CNR and spatial resolution). The digital system has the advantage of being easier to use than the system with films, with lower acquisition times. The use of a digital detector has been validated regarding the image quality necessary for characterization of bone mineralization. The next step concerns validation of quantitative assessment of bone mineralization variables as mean DMB, heterogeneity index. These variables will be measured and compared between the two systems for several bone specific units from the same human bone samples.

Fig. 1 .

 1 Fig. 1. Mass attenuation coefficients of aluminum (Al) and hydroxyapatite (HA) as a function of energy.

Fig. 2 .

 2 Fig. 2. Microradiographic systems using (a) films and (b) digital detector.

Fig. 5 .

 5 Fig. 5. Contrast to noise ratio as a function of high voltage either with a constant power (solid line) or constant exposure values (dotted line).

Fig. 6 .

 6 Fig. 6. Modulation transfer function curves obtained for several magnifications M with digital detector (full lines) compared to film (dashed line).

Fig. 7 .

 7 Fig. 7. Calculation of the contrast transfer functions (CTF). (a) image of the bar space pattern, (b) the corresponding grey-level profile and (c) CTF curves obtained for film and digital detector with different magnifications.

Fig. 8 .

 8 Fig. 8. (a) Comparison of microradiographs of a human iliac bone biopsy and zoom of a specific part of the image obtained with (b) the radiographic film and with the digital system at two different magnifications: (c) 3.7 and (d) 10.8. (e) Curves representing the grey-level profiles obtained respectively on the images (b), (c) and (d) across two bone specific units.

Table 1 .

 1 Technical specifications of both systems using radiographic films and a digital detector.

	Exposure	High voltage: 25 kV	High voltage: 40 kV
	parameters	Current: 25 mA	Current: 50 µA
		Power: 625 W	Power: 2 W
	Source to object	About 30 cm	Range: 1-25 cm
	distance
	Object to detector	Contact	Range: 1-25 cm
	distance
	Detector	Kodak radiographic films	Photonic science CCD camera FDI VHR 11M
		resolution: over 1000 pl/mm	active area: 36x24 mm (4008 x 2671 pixels)
		Revelator: Kodak D19	scintillator: Gd2O2S:Tb
		Fixator: Ilford Hypam	filter Al: 12 µm
		Temperature: 20°C
	Image digitization	Film digitized using a microscope	12 bits digital image
	step	coupled with a digital camera with 8	Detector pixel size: 9 µm
		bits coding
			Object pixel size depends on magnification
			factor.

Features

System using radiographic films Digital detector system X-ray tube X ray generator: Philips PW 1830/40 Diffraction tube: Philips PW 2273/20 Power maxi: 2.2 kW Cu anode, Ni filter Beryllium window: 300 µm Focal spot size: 4.2 mm 2 Microfocus Hammamatsu X-ray system L9421 -02 Power: maxi 8 W Cu anode, Ni filter Beryllium window: 150 µm Focal spot size: 5 µm diameter

Table 2

 2 Acquisition settings for digital system

	Parameter	Value
	High voltage (kV)	40
	Current (µA)	50
	Acquisition time (s)	7
	Number of frames	6
	Magnification factor	10.8

Acknowledgements 13

This work was supported by the "Région Rhône-Alpes" (CIBLE N° 11 010877 01) and a grant "Bonus Qualité Industrielle: QUATIOS" OSEO ANVAR RHONE-ALPES 2010. Part of this work was conducted in the frame of the ANR project MALICE N° BS09-032, and within the framework of the labex CeLyA of Université de Lyon operated by the French National Research Agency (ANR-10-LABX-0060/ANR-11-IDEX-0007). The authors thank Perceval Guillou for software assistance during installation of the system and testing, David Marchat (Assistant Professor, Center for Biomedical and Healthcare Engineering, Ecole Nationale Supérieure des Mines de Saint-Etienne) for providing us a hydroxyapatite sample and Jean-Paul Roux (UMR1033) for technical and scientific discussions.


Bertrand Payet
email: bpayet@univ-reunion.fr

Alain Shum

Cheong Sing

Jacqueline Smadja

Potentialities of molasses from cane sugar crystallization in food formulations

The molasses gave "amber", "golden", "caramel-like" colour to the ice creams. Sweetness of the ice cream was considerably enhanced by the molasses and was very pleasant for the ice cream with A molasse. The sweet taste of the ice creams with B and C molasses was very intense and too excessive, formulas with smaller amounts of these products should be considered. A caramel flavour was given to the vanilla ice creams by the molasses and malaga note was often cited for the ice creams containing A and B molasses. The vanilla flavour of ice creams was masked by the use of molasse C. Vanilla ice cream with A molasse was the most appreciated for its colour and taste. ■ Vanilla ice cream with A molasse ▲ Vanilla ice cream with B molasse • Vanilla ice cream with C molasse ♦ Traditional vanilla ice cream

Introduction

Molasses are produced together with granulated sugar. Three kinds of molasses are obtained during fractional crystallization : A, B, C molasses issued from the first, the second and the third step of crystallization, respectively. Molasses are made up 40% to 60 % of sugars (mainly sucrose). [START_REF] Higgingotham | Sugar Technology[END_REF] Several parameters such as the sucrose inversion, the presence of aminoacids, water, and mineral salts, the alkaline pH, the sugar process duration (3 days) and the high temperature (90°C) favour Maillard reactions. These reactions generate coloured macromolecules and heterocycles which have a strong olfactive impact despite their low concentration and are responsible for the "grilled", "roasted", "caramel" and "burned" olfactive notes of molasses. [START_REF] Vernin | [END_REF] Vanilla ice creams, sweet buns, butter cakes were made according to standard formulas. 10% to 20% of sugar quantity were substituted by A, B, C molasses. The traditional formulas and the new formulas were compared by tests relating to physicochemical and sensorial parameters.

Conclusion

Physicochemical parameters (humidity loss, melting temperature kinetic of ice creams) were not significantly modified by the addition of molasses to the food.

Nevertheless, in all cases, molasses changed sensorial aspects of the prepared foodstuffs. Differentiation between the molasses was possible in the formulated foods. Among the molasses, the C molasse was the most enriched in Maillard reaction products having a strong impact on aromatic notes and colour. The A molasse had a fewer impact on aromatic notes and colour but led to preferred formulated foods.

Instead of their uses in animal foods or as substrate for alcoholic fermentations, molasses could be considered as interesting sources of flavourings.

LCSNSA

Sensory analysis applied to butter cakes Sensory analysis applied to sweet buns

Substitution of sugar by molasses in the butter cakes and in the sweet buns gave them flavour and better aromatic properties. The Maillard reaction products are responsible for the differences in aromatic profiles of these foodstuffs. Molasses gave candied and fruity notes to the baked products and also reminiscent of gingerbread.

Cakes and buns prepared from standard formulas were less roasted than those prepared with A, B and C molasses. The deep colour given by the B and C molasses was, however, moderately well appreciated. No differences were observed in term of mouthfeel for the four butter cakes. A and C molasses greatly improved the mouthfeel of the sweet buns. Thirty two compounds were identified in the headspace of the molasses including fourteen furans and pyrazines.

During the cane sugar process, an increase of furan and pyrazine contents was observed from A molasses to C molasses according to the evolution of sugar composition. 4

■■

  Butter cake with A molasse ▲ Butter cake with B molasse• Butter cake with C molasse ♦ Traditional butter cake Sweet bun with A molasse ♦ Traditional sweet bun• Sweet bun with C molasse ▲ Sweet bun with B molasse

31.6 111.0 95.9

	2,5-dimethylfuran		27.8 68.4
	2-butylfuran	8.9	43.7 46.0
	2-pentylfuran		79.8 72.1
	2-furancarboxaldehyde	19.7 59.9 79.8	caramel
	2-furanmethanol, acetate		4.8	42.8	smooth, floral
	1-(2-furanyl) ethanone			53.3
	2-furancarbothioic acid, S-methyl ester			80.2
	5-methylfurancarboxaldehyde		16.1 44.0	spicy, caramel
	2-furanmethanol	12.2 63.5 191.6	warm, burnt
	total 40.7 387.6 736.4
	Pyrazines and derivated
	2,5-dimethylpyrazine	12.7 37.9 24.7	peanut, roasted
	2,6-dimethylpyrazine	18.9 58.1 38.4	cocoa, roasted
	2-ethyl-6-methylpyrazine			18.7	peanut, roasted
	2-ethenyl-6-methylpyrazine		15.0 14.1
	total


Telesphore Tiendrebeogo

Damien Magoni
email: magoni@labri.fr

Virtual and Consistent Hyperbolic Tree: A New Structure for Distributed Database Management

This paper describes a new scalable, reliable and consistent structure for implementing a distributed database system. This structure is based on the hyperbolic geometry and we call it a Virtual and Consistent Hyperbolic tree (VCH-tree). This structure is intended for supporting queries over possibly large spatial datasets distributed over interconnected servers. The VCH-tree is comparable to the well-known R-tree structure, but it uses the hyperbolic geometry properties of the Poincaré disk model. It uses a distributed balanced Q-degree spatial tree that scales with data objects' insertions into potentially any number of storage servers through virtual hyperbolic coordinates. A user application manipulates the structure from a client node. The client can connect to the system through one of the database servers that is already in the VCH-tree. Messages are then routed towards the proper server by a greedy algorithm which is using the virtual hyperbolic coordinates attributed to each server. We have performed simulations to assess the efficiency and reliability of the VCH-tree. Results show that our VCHtree is consistent and has expected performances given the scalability and flexibility it provides to distributed database applications.

Introduction

We aim at indexing large data sets of spatial objects, each uniquely identified by an object identifier (OID) and stored in the VCH-tree with a given address. We define a scalable and reliable index that generalizes R-tree structure for distributed data structure [START_REF] Silberschatz | Database Systems Concepts[END_REF] that we call VCH-tree. A VCH-tree permits redundancy of object references, like R-tree [START_REF] Guttman | R-trees: A dynamic index structure for spatial searching[END_REF] or R*-tree [START_REF] Jansson | Constructing the r* consensus tree of two trees in subcubic time[END_REF]. The fundamental principle of our system is to map a large database object identifier space onto a set of servers in a deterministic and distributed way. Roughly, given an object key, the system is able to obtain the location of several database servers where are stored the corresponding values. In order to route the queries, each database server usually maintains the status of its connections to all other databases servers which increases drastically the number of messages exchanged, and this may constitute a severe scaling limitation. The same applies to the number of routing hops that must not grow too fast with the number of database servers in the system [START_REF] Idowu | Transactions-distributed database systems: Issues and challenges[END_REF]. Moreover, most of distributed database systems suffer from the lack of flexibility in storage queries (i.e., values' placement) involving consequently a heavy lookup traffic load on the lookup paths of the underlying nodes. We argue that we can address and overcome simultaneously all the aforementioned issues by maintaining a good trade-off between robustness, efficiency and system complexity. To this end, we promote an indexing system for distributed database relying on hyperbolic geometry [START_REF] Anderson | Hyperbolic geometry[END_REF]. In this paper we make the following contributions:

-We provide a new architecture for indexation in the distributed database system without any constraints. The database servers can connect arbitrarily to each other, the data objects can be inserted, updated or deleted without the cost of maintaining any global knowledge of the servers' topology. Our approach is based on the groundbreaking work of Kleinberg [START_REF] Kleinberg | Geographic routing using hyperbolic space[END_REF], that we have enhanced in order to manage a dynamic topology which is able to grow and shrink over time. -We propose a specific key based greedy routing system for storing the mapping of database object identifiers to the servers' addresses in the hyperbolic plane in a flexible and efficient way. Values are stored in order to avoid overloading a particular zone of the distributed system. Furthermore, storage and retrieving queries can be solved within O(logN ) hops. -To improve database object availability and access performance, our system which is akin to a Distributed Hash Table (DHT) system embeds a redundancy and caching mechanism that can be adjusted to obtain a good trade off between reliability and storage consumption. -We have carried out simulations to demonstrate the interest of our solution in terms of feasibility to build a distributed database system over a topology based on the hyperbolic plane.

The remainder of this paper is organized as follows. Section 2 gives a brief overview of the related previous work. Section 3 highlights some properties of the hyperbolic plane when represented by the Poincaré disk model. Section 4 defines the local addressing and greedy routing algorithms of our distributed VCH-tree system. Section 5 defines the binding algorithm of our VCH-tree system. Section 6 presents the results of our practicability evaluation obtained by simulations and we conclude in Section 7.

Related Work

Until recently, most of the spatial indexing design efforts have been devoted to centralized systems [START_REF] Gaede | Multidimensional access methods[END_REF] although, for non-spatial data, research devoted to an efficient distribution of large data sets is well-established [2] [3]. Many Scalable Distributed Data Structure (SDDS) schemes are hash-based, e.g., variants of LH* [START_REF] Jajodia | Lh*re: A scalable distributed data structure with recoverable encryption[END_REF], or use a Distributed Hash Table (DHT) [START_REF] Crainiceanu | Querying peer-to-peer networks using p-trees[END_REF] [START_REF] Lakshman | Cassandra: A decentralized structured storage system[END_REF]. Some SDDSs are range partitioned, from RP* based systems [START_REF] Karlsson | hqt*: A scalable distributed data structure for high-performance spatial accesses[END_REF] up to BATON [START_REF] Gaede | Multidimensional access methods[END_REF] most recently. There were also proposals for k-d partitioning, e.g., k-RP [START_REF] Litwin | k-rp*s: A scalable distributed data structure for highperformance multi-attribute access[END_REF] using distributed kd-trees for data points, or hQT* [START_REF] Karlsson | hqt*: A scalable distributed data structure for high-performance spatial accesses[END_REF] using quad-trees for the same purpose. Hambrusch and Khokhar [START_REF] Hambrusch | Maintaining spatial data sets in distributedmemory machines[END_REF] present a distributed data structure based on orthogonal bisection trees (2-d KD trees). Kriakov et al. [START_REF] Kriakov | Management of highly dynamic multidimensional data in a cluster of workstations[END_REF] describe an adaptive index method which offers dynamic load balancing of servers and distributed collaboration. The structure requires a coordinator which maintains the load of each server.

Hyperbolic Geometry

The model that we use in our system to represent the hyperbolic plane is called the Poincaré disk model. In the Poincaré disk model, the hyperbolic plane is represented by the open unit disk of radius 1 centered at the origin. In this specific model:

-Points are represented by points within this open unit disk.

-Lines are represented by arcs of circles intersecting the disk and meeting its boundaries at right angles.

In this model, we refer to points by using complex coordinates. An important property is that we can tile the hyperbolic plane with polygons of any sizes, called p-gons. Each tessellation is represented by a notation of the form {p, q} where each polygon has p sides with q of them at each vertex. There exists a hyperbolic tessellation {p, q} for every couple {p, q} obeying (p-2) * (q -2) > 4. In a tiling, p is the number of sides of the polygons of the primal (the black edges and green vertices in figure 1) and q is the number of sides of the polygons of the dual (the red triangles in figure 1). Our purpose is to partition the plane and address each node uniquely. We set p to infinity, thus transforming the primal into a regular tree of degree q. The dual is then tessellated with an infinite number of q-gons. This particular tiling splits the hyperbolic plane in distinct spaces and constructs an embedded tree that we use to assign unique addresses to the nodes. An example of such a hyperbolic tree with q = 3 is shown in figure 1.

In the Poincaré disk model, the distances between any two points z and w are given by curves minimizing the distance between these two points and are called geodesics of the hyperbolic plane. To compute the length of a geodesic between two points z and w and thus obtain their hyperbolic distance d H , we use the Poincaré metric which is an isometric invariant given by the formula:

d H (z, w) = argcosh(1 + 2 × |z -w| 2 (1 -|z| 2 )(1 -|w| 2 ) ) ( 1
)
This formula is used by the greedy routing algorithm shown in the next section. We now explain in this section how we create the hyperbolic addressing tree for database servers interconnections and how queries can be routed in our distributed database system. The first step in the creation of a VCH-tree of servers nodes is to start the first database server and to choose the degree of the addressing tree. We recall that the hyperbolic coordinates (i.e., a complex number) of a server node of the addressing tree are used as the address of the corresponding database server in the distributed data base system. A server node of the tree can give the addresses corresponding to its children in the VCH-tree. The degree determines how many addresses each database server will be able to give for news nodes servers connections. The degree of the VCH-tree is defined at the beginning for all the lifetime of the distributed database system. The distributed database system is then built incrementally, with each new data server joining one or more existing data servers. Over time, the data servers will leave the overlay until there is no server left which is the end of the distributed database system. So, for every data object that must be stored in the system, an OID is associated with him and map then in key-value pair. The key will allow to determine in which data servers the object will be stored (as explained in the following section). Furthermore when a data object is deleted, the system must be able to update this operation in all the system by forwarding query. This method is scalable because unlike Kleinberg [START_REF] Kleinberg | Geographic routing using hyperbolic space[END_REF], we do not have to make a two-pass algorithm over the whole distributed system to find its highest degree. Also in our system, a server can connect to any other server at any time in order to obtain an address.

The first step is thus to define the degree of the tree because it allows building the dual, namely the regular q -gon. We nail the root of the tree at the origin of the primal and we begin the tiling at the origin of the disk in function of q. Each splitting of the space in order to create disjoint subspaces is ensured once the half spaces are tangent ; hence the primal is an infinite q-regular tree. We use the theoretical infinite q-regular tree to construct the greedy embedding of our q-regular tree. So, the regular degree of the tree is the number of sides of the polygon used to build the dual (see figure 1). In other words, the space is allocated for q child database servers. Each database server repeats the computation for its own half space. In half space, the space is again allocated for q -1 children. Each child can distribute its addresses in its half space. Algorithm 1 shows how to compute the addresses that can be given to the children of a database server. The first database server takes the hyperbolic address (0;0) and is the root of the tree. The root can assign q addresses. Algorithm 1 Calculating the Coordinates of a Server's Children.

1: procedure CalcChildrenCoords(server, q) 2:

step ← argcosh(1/sin(π/q)) 3:

angle ← 2π/q 4: childCoords ← server.Coords 5:

for i ← 1, q do 6:

ChildCoords.rotationLef t(angle) 7:

ChildCoords.translation(step) 8:

ChildCoords.rotationRight(π) 9:

if ChildCoords = server.P arentCoords then 10:

StoreChildCoords(ChildCoords) 11:

end if 12:

end for 13: end procedure This distributed algorithm ensures that the database servers are contained in distinct spaces and have unique coordinates. All the steps of the presented algorithm are suitable for distributed and asynchronous computation. This algorithm allows the assignment of addresses as coordinates in dynamic topologies. As the global knowledge of the distributed database system is not necessary, a new server can obtain coordinates simply by asking an existing server to be its parent and to give it an address for itself. If the asked server has already given all its addresses, the new server must ask an address to another existing database server. When a new server obtains an address, it computes the addresses (i.e., hyperbolic coordinates) of its future children(The new database servers which shall connect to the distributed database system). The addressing VCH-tree is thus incrementally built at the same time than the distributed database system.

When a new database server has connected to database servers already inside the distributed database system and has obtained an address from one of those database servers, it can start sending requests to store or lookup database object in the distributed database system. The routing process is done on each database server on the path (starting from the sender) by using Algorithm 2, a greedy algorithm based on the hyperbolic distances between the servers. When a query is received by a database server, the database server computes the distance from each of its neighbors to the destination and forwards the query to its neighbor which is the closest to the destination (destination database server computing is given in Section 5).

Algorithm 2 Routing a Query in the Distributed Database System. In a real network environment, link and server failures are expected to happen often. If the addressing VCH-tree is broken by the failure of a database server or link, we flush the addresses attributed to the servers beyond the failed server or link and reassign new addresses to those servers (some servers may have first to reconnect to other servers in order to restore connectivity). But this solution is not developed in this paper.

In this section we explain how our distributed database system computes the destination database servers addresses for storing and retrieving queries. Indeed, the first server contacted by a client (prime server) for sending a query in the system consider the latter as a data object that can be stored or looked up. Thus this server generates an OID associated to the data object and the latter is mapped onto hyperbolic addresses corresponding to destination database servers' addresses in the VCH-tree.

On startup, each new client query is associated with the data object with OID corresponding to the name of the query and that identifies the query it runs on. This name will be kept by data object during all the lifetime of the distributed database system.

When the prime database server computes some specific addresses of database servers, when it is about a storage query, it stores the name (OID) and value of query in these specific addresses of distributed database servers, thus the data object in the DHT, when it is about a retrieving query, it contacts database servers which addresses has been computed. In our distributed system, the name is used as a key by a mathematical transformation. If the same name is already stored in the distributed database system, an error message is sent back to the prime server (Server by whom the client is directly bound) in order to generate another name. Thus the distributed database system structure itself ensures that OIDs are unique.

A (OID, value) pair, with the name acting as a key by mapping is called a binding. Figure 2 shows how and where a given binding is stored in the distributed database system. A binder is any database server that stores these pairs. The depth of a server in the addressing VCH-tree is defined as the number of parent servers to go through for reaching the root of the VCH-tree (including the root itself). When the distributed database system is created, a maximum depth for the potential binders is chosen. This value is defined as the binding VCH-tree depth. To ensure a load balancing of the system, the depth p is chosen such p minimize the inequality 2 :

p × ( (1 -(p -1) d ) 2 -p + 1 ≥ N (2)
with p the depth, d the degree and N the number of database servers.

When a new database server joins the distributed database system by connecting to other servers, it obtains an address from one of these servers. Next, the server stores its own binding in the system. So, during his life, each database server tries to join others by sending a join query. Each server cannot accept that a limited number of join queries independently of the degree of the VCHtree. The new connections serve as shortcuts during the phases of storage and retrieving of data objects. We call these connections, shortcut links as indicated in Figure 2.

Storage Query Processing

When a client wants to send a storage query (i.e., insertion), the first server with whom it is connected consider a query as an object (thus generating an OID) and creates a key by hashing its name with the SHA-512 algorithm. It divides the 512-bit key into 16 equally sized 32-bit subkeys (for redundant storage). The server selects the first subkey and maps it to an angle by a linear transformation. The angle is given by:

α = 2π × 32-bit subkey 0xFFFFFFFF (3)
The database server then computes a virtual point v on the unit circle by using this angle:

v(x, y) with x = cos(α) y = sin(α) (4)
Next the database server determines the coordinates of the closest binder to the computed virtual point above by using the given binding tree depth.

In the figure we set the binding VCH-tree depth to three to avoid cluttering the figure. It's important to note that this closest binder may not really exist if no database server is currently owning this address. The database server then sends a storage query to this closest database server. This query is routed inside the distributed database system by using the greedy algorithm of Section 4. If the query fails because the binder does not exist or because of database server/link failures, it is redirected to the next closest binder which is the father of the computed binder.

The path from the computed closest binder to the farthest binder is defined as the binding radius. This process ensures that the queries are always stored first in the binders closer to the unit circle and last in the binders closer to the disk center. However to avoid overloading the farthest binder particularly and to keep a load balancing, we limit the number of storages S like follow:

S ≤ 1 2 × log(N ) log(q) (5)
with N equal to number of distributed database servers, q to degree of VCHtree.

Furthermore the previous solution, any binder will be able to set a maximum number of stored queries and any new database server to store will be refused and the query redirected as above. Besides, to provide redundancy and so ensure the availability and reduce the latency period in the lookup process, the database server does the storage process described above for each of the other 15 sub keys. Thus 16 different binding radii will be used at the most and this will improve the even distribution of the pairs (key-value).

In addition to this and still for redundancy purposes, a pair key-value of the data object may be stored in more than one database server of the binding radius. A binder could store a data object and still redirect its query for storage it in other ancestor binders. The number of stored copies of a key-value pair along the binding radius may be an arbitrary value set at the distributed system creation. Similarly the division of the key in 16 sub keys is arbitrary and could be increased or reduced depending on the redundancy needed. To conclude we can define two redundancy mechanisms for storage copies of a given binding:

1. We can use more than one binding radius by creating several uniformly distributed subkeys. 2. We can store the data object key-value pair in more than one binder in the same binding radius.

These mechanisms enable our distributed database system to cope with a non-uniform growth of the database servers and they ensure that a data object will be stored in a redundant way that will maximize the success rate of its retrieval. The numbers of sub keys and the numbers of copies in a radius are parameters that can be set at the creation of the distributed database system. Increasing them leads to a tradeoff between improved reliability and lost storage space in binders. Besides our solution has the property of consistent hashing: if one database server fails, only its keys are lost but the other binders are not impacted and the whole system remains coherent. Algorithm 3 illustrates the previous mechanism. OID ← Query.GetOID() 3:

Key ← Hash(OID) 4:

for all red ∈ R Circular do 5:

depth ← PMax 6:

i ← 1 7:

while i ≤ 1 2 × log(N ) log(q)
&& depth ≥ 0 do 8:

SubKey[red][depth] ← ComputeSubkey(Key)[red][depth] 9: T argetServerAddr[red][depth] ← ComputeAddr(SubKey[red][depth]) 10: T argetServer ← GetT arget(T argetServerAddr[red][depth]) 11:
if route(Query, T argetServer) then 12:

i + + 13:

put(OID, Query) 14:

end if 15:

depth --16:

end while 17:

end for 18: end function

Lookup Query Processing

Now, if the client wants to lookup a data object in the distributed database system, a prime server is contacted and generates an OID for the client query. Here again, the OID is mapped into a key by the SHA-512 algorithm, thus the 512 bits key is divided into 16 subkeys. Each subkey, by the process described in Section 5.1, will be transformed into an address that represents the address of the database server where the data object is stored. The latter is contacted by the prime database server for updating, deleting or retrieving the associated value. When the redundancy mechanism has been used to store the data object, the lookup repeats the latter process of lookup for any subkey, thus the operation will be performed on all database servers that contain the data object. Our distributed system ensures the coherence of data objects of the distributed database. This mechanism is illustrated by Algorithm 4.

Evaluation

We performed experiments for evaluating the behavior of a VCH-tree over large datasets. Furthermore, we consider that the system is static so there are no Algorithm 4 Lookup and Update Algorithm.

1: function Lookup(Query) return V alue 2:

QueryOID ← T arget.GetQueryOID() 3:

Key ← Hash(QueryOID) 4:

for all red ∈ R Circular do 5:

depth ← PMax 6:

i ← 1 7:

while i ≤ 1 2 × log(N ) log(q)
&& depth ≥ 0 do end for 25: end function join or leave of database servers during the simulation. We use the Peersim [START_REF] Montresor | PeerSim: A scalable P2P simulator[END_REF] simulator for running event-driven simulations. The study involved the following parameters of the VCH-tree:

-The number of database servers connected and used to store the data objects.

Here we have considered 10000 database servers; -Each run lasts 2 hours of simulated time; -We try to store 6 millions of data objects in our distributed database system following an exponential distribution with a median equal to 10 minutes; -The maximum capacity for each server is set to 6000 objects.

We studied the behavior of our structure for both data objects' storage and retrieval in the system. We are interested in observing the scalability of our system, the shape of the hyperbolic tree, the storage load balancing and the length of the paths of the queries.

Spatial Shape of the VCH-tree

Figure 3 shows an experimental distribution of points corresponding to the scatter plot of the distribution of the database servers in our system. We can see that our VCH-tree is balanced. Indeed, we can notice by part and others around the unit circle which we have database servers. This has an almost uniform distribution around the root, which implies that our system builds a well-balanced tree that will more easily allow to reach a proper load balancing for the storage. Figure 4 shows correspondingly Poincaré disk model that no address of database server belongs on the edge of the unit circle. Indeed, the addresses of database server were obtained by projection of the tree of the hyperbolic plane in a circle of the Euclidian plane of radius 1 and of center with coordinates (0; 0). This result shows that our distributed database system can grow towards infinity in theory. In practice, other parameters such as real number precision do bring limitations.

Load Balancing in the VCH-tree

Figure 5 shows a plot of the average number of objects stored by the database servers over time. So this figure shows a regular growth of this number of data objects stored in function of time. Indeed, 293.27 data objects on average are stored by database server after 10 minutes vs 620.4 after 2 hours. It is interesting to notice that the standard deviation remains low, approximately at 10 % of the average. This indicates a low dispersal of the number of objects stored on the servers during the simulation.

Indeed, if we use our results to build the confidence interval, we can say that after 10 minutes of simulation, 68.2 % of the database servers store between 263.69 and 322.71 data objects and 95 % store between 234.18 and 352.22 against 68.2 % of the database servers which store between 560.18 and 681.58 data objects after 2 hours and 95 % of the database servers who store between 497.95 and 742.84 data objects after 2 hours. In view of these results, we can say that our system maintains a proper load balancing between database servers which ensures the stability of our distributed database system.

Storage and Retrieval Efficiency in the VCH-tree

Figure 6 and 7 show that during the simulation, queries in both cases can be answered within O(log N) where N is equal to the number of database servers in the system. As the standard deviation is very low (less than 5 % of the average for storage and retrieval), we did not represent it on the figures. In the worst case, queries need to travel less than 4 database servers in the system for either storing, or retrieving a data object. Besides, what is also interesting to note is that the plot decreases slowly to become stationary after around 100 minutes in both cases. It can be explained because during the simulation, the database servers create shortcuts as indicated in Section 5. These shortcuts allow to reach their target in fewer hops. The stationary situation is understandable by the fact that after a while, all the database servers reached their maximum number of shortcuts created and the most part of the queries is processed on average in less than 3.75 hops in both cases.

Conclusion

In this paper, we have presented a new structure called VCH-tree. This hyperbolic tree presents some properties that allow us to propose a consistent system of distributed database servers using virtual addresses made from hyperbolic coordinates. We have evaluated the performances of our system by simulation. We have shown that our system is scalable in terms of the number of database servers that can be interconnected as well as in terms of the number of hops to route the queries. We have also shown that the placement of the different database servers allows us to keep a well-balanced tree. Furthermore, we have shown that our system maintains a load balancing for the storage of data objects. For future work, we plan to study our solution comparatively to the other ones described in the state of the art in order to assess its benefits relatively to those existing solutions.

Fig. 1 .

 1 Fig. 1. 3-regular tree in the hyperbolic plane

Fig. 2 .

 2 Fig. 2. Storage in the VCH-tree

Algorithm 3

 3 Storage Algorithm. 1: function Store(Query) 2:

Fig. 3 .

 3 Fig. 3. Scatter plot of the spatial positions of the database servers.

Fig. 4 .

 4 Fig. 4. Scatter plot of the positions of the servers in the neighborhood of the unit circle.

Fig. 5 .

 5 Fig. 5. Average load on the database servers over time

Fig. 6 .Fig. 7 .

 67 Fig. 6. Path length of storage queries