G Sutter

S Philippon

F Garcin

Dynamic analysis of the interaction between an abradable material and a titanium alloy

Keywords: Abradable material, Blade/casing interaction, High-speed experimental device, Force measurement

The purpose of the present paper is to study interaction phenomena between an abradable material (M601) and a titanium alloy (TA6V). This interaction which can occur in an aircraft engine when the rotating blade tips rub on the compressor casing, is analysed on an experimental laboratory set-up. Only orthogonal cutting conditions are considered in this preliminary study to monitor the load evolution with velocity. The experimental results obtained for a range of velocities 10 < V < 107 m/s and for a range interaction depth 0.05 < t 1 < 0.5 mm are presented. After analysis of measurements, the mean interaction force increases proportionally to the interaction depth increment. The effects of the contact geometry on the loads, especially the radius of the edge, are investigated. The increase of the radius value for a given interaction depth increases the load value.

This work provides new basic data for behaviour laws introduced in numerical models simulating the interaction between the blade tip and the compressor casing in turbine engine applications. Original high-speed photographs recorded in real time during the process illustrate the chip separation by brittle manner. A temperature field in the abradable material will be able to be obtained for only a low speed of 5 m/s.

Introduction

The optimisation of compressor engines is of prime importance for the improvement of aircraft performances. The increasing compressors efficiency and therefore the reduced fuel consumption can be obtained by minimizing the clearance between compressor casing and the rotating blade tip [START_REF]Abradable seal with improved properties[END_REF][START_REF] Dorfman | Gas turbines use abradable coatings for clearance-control seals[END_REF]. The adopted technical solution consists in the use of the blade tip, as a cutting tool, to scrape the casing. This method is widespread in other turbine-engine applications such as labyrinth seals in bearings [START_REF] Mutasim | Evaluation of plasma sprayed abradable coatings[END_REF] or the turbocharger of diesel engines [START_REF] Sato | Development of the new TD27Eti engine for recreational vehicles[END_REF].

Therefore the compressor casing is covered with abradable material. This coating, composed of metal phase and selflubricating non-metal phase with high porosity, offers a good balance between abradability and erosion resistance [START_REF] Bounazef | The wear, deterioration and transformation phenomena of abradable coating BN-SiAl-bounding organic element, caused by the friction between the blades and the turbine casing[END_REF][START_REF] Maozhong | Erosion wear behaviour and model of abradable seal coating[END_REF][START_REF] Maozhong | Friction and wear behaviour and abradability of abradable seal coating[END_REF]. In addition, the smooth surface obtained after the interaction pro-cess, favourable to aerodynamics properties, must be generated without damage to the blade tip.

Unfortunately, turbine rotor misalignment, start-ups, stops, vibrations, thermal and centrifugal effects on the rotating blades create undesirable interaction. The casing geometry is then modified, increasing the clearance and damaging the blades. Particles of abradable material of the coating are transferred onto the blades increasing the mass and are ingested by the engine [START_REF]Abradable seal with improved properties[END_REF]. The necessity for the aircraft engineering to understand the interaction behaviour allows to avoid premature turbine wear or failure, and then to increase the useful lives of engines.

The wear mechanisms occurring in abradable seals of gas turbines have been identified by an investigation of worn components, see Borel et al. [START_REF] Borel | The wear mechanisms occurring in abradable seals of gas turbines[END_REF]. The main interaction mechanisms are cutting, deformation, adhesive transfer and melting. Three other mechanisms occur without contact: erosion, corrosion and heating. In these last cases, roughness, surface chemistry and color have been changed. The deformation process creates a striated surface and compacts coating leading to more future violent interactions. The adhesive transfer corresponds to the transfer of material from the coating onto the blade or reciprocally and is often associated with the deformation mechanism. This material transfer may form an alloy having a lower melting point than the original constituents. The melting of these alloys produces a hard layer during the following rubbing interaction and can seriously damage the blade. Thus, deformation, adhesive transfer and melting are not favourable and must be avoided. When the cutting interaction is observed, a new smooth surface with a shiny appearance is generated and the blade tip is not damaged.

Extensive research about erosion resistance and abradability has been carried out [START_REF] Bounazef | The wear, deterioration and transformation phenomena of abradable coating BN-SiAl-bounding organic element, caused by the friction between the blades and the turbine casing[END_REF][START_REF] Maozhong | Erosion wear behaviour and model of abradable seal coating[END_REF][START_REF] Maozhong | Friction and wear behaviour and abradability of abradable seal coating[END_REF]. The abradability tests were carried out on a rub test machine [START_REF] Maozhong | Friction and wear behaviour and abradability of abradable seal coating[END_REF] for small sliding velocity V = 0.377 m/s and for loads varying from 40 up to 100 N. The results obtained show the abradability decrease with the increase of hardness for a given kind of coating (the abradable specimens tested were M307, M310, M313 and M601). But for similar hardness of different kinds of coating, abradability is very different. So, the choice of the coating cannot be confirmed only with hardness. The erosion wear behaviour was investigated by a sand erosion machine for impact velocity in a range from 30 to 80 m/s and for an impact angle in a range from 30 • to 90 • [START_REF] Maozhong | Erosion wear behaviour and model of abradable seal coating[END_REF]. The erosion rate is related to the characteristics of the abradable material such as composition, microstructure and hardness, but it is not related to the thickness of the coating. At 30 • impact angle, micro-cutting plowing and tunneling of the non-metal phase was observed.

But the frequently used test device is the Sulzer Metco rig. It consists of a rotor on which a blade is mounted while the abradable material specimen is fixed on a dynamometer. The temperature measurements during the process can be obtained by a thermocouple inserted in the specimen [START_REF] Bounazef | The wear, deterioration and transformation phenomena of abradable coating BN-SiAl-bounding organic element, caused by the friction between the blades and the turbine casing[END_REF] or a thermocouple array fixed on the back surface of the blade [START_REF] Laverty | Rub energetic of compressor blade tip seals[END_REF] or by optical pyrometers [START_REF] Laverty | Rub energetic of compressor blade tip seals[END_REF][START_REF] Wang | Criteria for analysis of abradable coatings[END_REF][START_REF] Emery | An experimental investigation of temperatures due to rubbing at the blade-seal interface in an aircraft compressor[END_REF]. The effects on the interaction mechanism of the blade tip incursion rate (from 2 to 2000 m/s), incursion depth (a few tenths of millimeter) blade tip velocity (from 100 to 450 m/s), coating temperature, coating properties, rubbing duration, blade material and shape can be investigated [START_REF] Dorfman | Gas turbines use abradable coatings for clearance-control seals[END_REF]. Note that generally, the edge of the blade tip is at right angle with motion direction and parallel to the abradable specimen [START_REF] Emery | An experimental investigation of temperatures due to rubbing at the blade-seal interface in an aircraft compressor[END_REF]. This rig has been used to analyse the transfer of coating onto the blade for a blade velocity in a range from 200 to 400 m/s and for incursion velocity in a range from 10 to 500 m/s [START_REF] Bounazef | The wear, deterioration and transformation phenomena of abradable coating BN-SiAl-bounding organic element, caused by the friction between the blades and the turbine casing[END_REF]. The obtained results show that if the blade incursion is set deep into the abradable coating, the material transfer is lessened. In the same way, an increase of the blade velocity decreases the coating material transfer onto the blade. The normal and transverse forces due to the interaction between different coatings sprayed on the blade tip and seal material (80% Ni-20% Cr) were measured with the same rub device where the blade tip velocity was maintained constant at 100 m/s [START_REF] Emery | An experimental investigation of temperatures due to rubbing at the blade-seal interface in an aircraft compressor[END_REF]. For incursion rates of 0.0127 and 0.0254 mm/s, with an incursion depth of 0.5 mm, the normal and transverse force varied approximately from 200 to 520 and 50 to 150 N, respectively. In addition the temperature measurements of the abradable seal and the blade were obtained during the interaction process by two optical pyrometers [START_REF] Emery | An experimental investigation of temperatures due to rubbing at the blade-seal interface in an aircraft compressor[END_REF]. In these tests, the transverse force and temperature, which varied from 500 to 1170 • C, seemed to be reasonably well correlated. This was confirmed as this temperature was reached due to the frictional energy which is directly proportional to the transverse force. To investigate temperature at an incursion rate of 0.75 mm/s, two optical pyrometers were used during the interaction process of four titanium blades rubbing at 152.4 m/s against two kinds of sealling material (NiCrAl-silicate and AlSi) [START_REF] Wang | Criteria for analysis of abradable coatings[END_REF].

In the case of titanium blade and NiCrAl-silicate seal interaction, the measured temperatures were 1259 and 881 • C, respectively. These temperatures were lower (601 and 798 • C, respectively) when the titanium blade was rubbed against the AlSi seal.

A different configuration was used to determine the interaction temperature with the following parameters: incursion rate from 0.0025 to 0.25 m/s and blade tip velocity from 152 to 213 m/s [START_REF] Laverty | Rub energetic of compressor blade tip seals[END_REF]. In these tests, the abradable material (hastelloy × feltmetal fibermetal) was applied on the rim of a rotating disk and the instrumented blade was mounted on the dynamometer. The blades made of titanium or Incoloy 901 were skewed at 45 • from the plane of the disk rotation to reproduce rub geometry similar to the blade tip of the compressor. The maximal temperature reached was about 1300 • C, although a quick quenching follows this increase. The incursion rate was found to be the most important independent variable while the blade tip velocity and the blade thickness showed moderate influence.

To sum up, the coating/blade interaction and results are important in the engine efficiency. This complex process is influenced by many parameters and gets complicated with the need for minimal clearances combined with the high speed. The extreme friction conditions between abradable coating and titanium blade can become identified with the domain of high-speed machining. In this work, an experimental approach of the interaction mechanism is developed by means of an original device used for the high-speed machining study [START_REF] Sutter | An experimental study of high-speed orthogonal cutting[END_REF]. The particular behaviour of the abradable material during the cut (weak cutting force), required the development of a sensitive load sensor. Two different solutions were compared in the presented results to define the most adapted to the interaction conditions. The influence of the interaction depth is analysed for a range of interaction speed reaching 110 m/s. Directly connected to manufacture, the edge radius of the blade tip is studied for two values corresponding to rectification and to taking the burrs off. To complete this work, high-speed photographs were taken in real time during the interaction process. A temperature field in the abradable material will be able to be obtained for only a low speed of 5 m/s.

Experimental device

The interaction between a compressor blade and the casing is analysed by means of an experimental device using a pneumatic gun. This propulsion device has been designed to obtain high-speed orthogonal cutting behaviours [START_REF] Sutter | An experimental study of high-speed orthogonal cutting[END_REF]. The maximal speed obtained with this device can reach 120 m/s. Although the velocities reached at the compressor blade extremities are about 500 m/s for the first compressor stage. This paper presents the work of a pre-study allowing the validation of the measurement chains at high speeds. Of course the analysis perspectives are to expand the explored velocities range in order to reach the real velocities. This experimental device is mainly composed (see Fig. 1) of a pneumatic air gun using the compressed gas expansion to launch a projectile. This projectile supports workpiece A which is covered with abradable material, and is guided along a launch tube during its acceleration step. The accurate tolerances of the tube geometry and the good surface condition allow to obtain an optimal precision of the projectile trajectory and ensure the workpiece position during the interaction with specimen B. Projectiles with different densities are used in terms of the velocity ranges and forces encountered in order to ensure quasi-stationary conditions. Note that the kinetic energy decrease during the interaction corresponds to a maximal relative velocity variation of 4%. The high speed (about 120 m/s) is obtained with a projectile made of nylon.

A second tube (receiving tube) allows to guide the projectile towards a shock absorber designed to receive the workpiece without damage for post-experiment analysis. The sideextremities of the parallelepiped workpiece A are covered with abradable material whose chemical composition is presented in Table 1. This type of material is only characterized by the Rockwell hardness (HR15y) for which the value is 51 [START_REF] Maozhong | Friction and wear behaviour and abradability of abradable seal coating[END_REF]. Specimen B made of a titanium alloy (Table 2) simulates the behaviour of the engine rotor blades. They are symmetrically fixed on a load sensor. The interaction depth between the abradable material and specimen B is modulated using the combination of the distance separating specimen B and the length of workpiece A. A groove is machined into the support (Fig. 2) before the abradable material deposit in order to minimize its unsticking during the interaction. Fig. 2 presents a specimen A with only one extremity covered with abradable material to show the preparation before the material deposit.

A set of photo-diodes placed at the launching tube exit and coupled with time counters establishes the record of the velocity.

The load sensor, detailed in Fig. 3, supporting specimen B is equipped with strain gauges in order to measure the axial forces generated by the contact between the abradable and the titanium alloy. This sensor has been calibrated statically and dynamically by means of impacts at low speed. A double validation of the strains measured by this sensor has been performed. Firstly, Fig. 2. Specimen A. numerical analysis by finite elements confirmed the calibration. Secondly, a KISTLER piezoelectric strain sensor (piezoelectric sensor for indirect measurement) was implanted in this support to check the set sensitivity. The evolution of sensors exit signals presented in Fig. 4 confirms important sensitivity of both load sensors and in addition their linearity in terms of load.

Different tests were carried out at high interaction velocities using a first load sensor equipped with strain gauges and a second with a piezoelectric sensor. The signals, obtained respectively by the strain gauges and the piezoelectric sensor, are presented in Figs. 5 and6. Although the sensitivity of these two sensors remains close, the frequency of the piezoelectric sensor due to the mass disrupts the measurement by superimposition of a damped oscillation. These disturbances are all the more intense as the interaction speed is high.

In spite of damped oscillations visible on the recording obtained with the piezoelectric sensor, the mean values of the two signals correspond to the same interaction forces. These tests confirmed the calibration of the sensor using deformation gauges. The results presented below are obtained with the second sensor.

To analyse the behaviour of the two materials during the process, a high-speed camera is located near the zone of interaction. Due to the high velocities of the interaction, the phenomenon duration time to be filmed requires very short exposure time (of the order of a few microseconds). Significant luminous intensity is also needed, which is ensured by two flashes of high power. This video recording requires a perfect synchronization of the trigger mechanism with the photographic recording system. Fig. 5. Typical signal obtained with strain gauges located on the load sensor. Depth of interaction t 1 = 0.10 mm, interaction velocity V = 44 m/s, edge sharpness r = 0.001 mm. Fig. 6. Typical signal obtained with KISTLER piezoelectric sensor located in the load sensor. Depth of interaction t 1 = 0.10 mm, interaction velocity V = 37.2 m/s, edge sharpness r = 0.001 mm.

Results and discussion

The first results presented in this paper are only a simplified approach of the real phenomena. Indeed an orthogonal interaction is simulated with our device in order to understand the intricate mechanisms which are taking place in the blades/casing interface. The edge of the moving specimen B is perpendicular to the translation direction. It is possible to establish a similarity with an orthogonal cutting process (Fig. 7). Furthermore the inclined plane of the profile of specimen B is 6 • comparable to a clearance angle. Various edge sharpness, noted r in Fig. 7, are imposed to specimen B to study the influence of the finishing process during the compressor blade manufacture. The tests pre- sented here have been performed with two edge radii r = 0.001 and 0.03 mm.

Fig. 8 presents the evolution of the measured forces F, created during the interaction between the two materials M A (abradable material) and M B (titanium alloy) for a translation velocity reaching 107 m/s and an edge radius r = 0.001 mm. An increase in forces is observed with the increase of the translation velocity. This sensitivity is at its highest when the interaction depth is high.

The intensity of the force does not seem to be linearly following the interaction depth t 1 . Indeed, for a velocity V = 50 m/s and t 1 = 0.05 mm, the average force F is 65 N. At the same velocity, an interaction depth t 1 = 0.1 mm corresponds to a force of 112 N and for t 1 = 0.5 mm, F = 405 N. However we have to assume that the contact localized under the edge of specimen B generates a residual force. This residual force generated by friction, exists even with an interaction depth t 1 = 0. This phenomenon has been experimentally verified for a decrease of the interaction depth close to 0, the interaction force reaches a minimum value different from 0. This residual force, noted F r , depends on the interaction speed but is mainly influenced by the value of the edge radius. Following this analysis, an edge radius of 0.001 mm combined with a minimal depth t 1 leads to a force F r = 50 N for the lower velocities and 65 N for higher velocities.

For some edge radius and some interaction speed, the measured force F can be also defined by the amount of a variable force F i and a constant force F r (see Fig. 7). This variable inter- action force can be represented by the following equation:

F i = F -F r ( 1
)
Thus for t 1 = 0.1 mm, the interaction force F i = 112 -50 = 62 N and for t 1 = 0.5 mm, F i = 355 N. These forces F i are proportional to the interaction depth.

In order to study the influence of the edge radius on the evolution of the forces, the previous tests are performed with a different value of r = 0.03 mm. The increase of the radius leads to the increase of interaction forces as seen in Fig. 9. The increasing force with respect to the edge radius is due mainly to the increase of the friction forces at the interface between two materials. Indeed a higher radius increases the friction area which leads to decrease the contact pressure. Consequently the friction coefficient increases as well as the residual force F r .

The interaction area is detailed in Fig. 10 in static position and Fig. 11 presents the video recording during the interaction process. As expected in the jet engine, the interaction between the titanium alloy and abradable material does not generate con- tinuous chips or too big particles. These pictures confirm that the interaction produces fine particles obtained by brittle manner already observed by Borel et al. [START_REF] Borel | The wear mechanisms occurring in abradable seals of gas turbines[END_REF].

An infrared video camera has been used to record the temperature gradient during the interaction. Because of the low dynamic of this video camera it has not been possible to record at high velocity. On the other hand we have used the high-speed numerical camera as a pyrometer to detect the temperature field near the interaction zone. This original technique was developed by Sutter et al. [START_REF] Sutter | An experimental technique for the measurement of temperature fields for the orthogonal cutting in high speed machining[END_REF] and Ranc et al. [START_REF] Ranc | Temperature measurement by visible pyrometry-orthogonal cutting application[END_REF] to record the temperature during high-speed machining. Unfortunately, the evolution of the temperature during the interaction is not large enough to be measured by this set-up. In spite of the very good resolution of this measurement device, the hot zones may be too localized. In addition the emission of particles at high speeds may obstruct the camera field of vision.

Fig. 12 presents the temperatures reached with a velocity of about 5 m/s. The low increase of the observed temperature predicts a limited influence of the temperature during the interaction under these conditions.

The reduced interaction velocity imposed by limited frame rate (200 Hz) of the infrared camera can explain the low measured temperature. Rapid increasing temperatures with the velocity allow to reach similar high temperatures as those previously mentioned [START_REF] Wang | Criteria for analysis of abradable coatings[END_REF][START_REF] Emery | An experimental investigation of temperatures due to rubbing at the blade-seal interface in an aircraft compressor[END_REF]. Moreover, in the previous works, the experiments are carried out on a rub test machine where the interaction process is repeated several times during the test. Consequently these repetitive contacts extend the thermal affected zone by conduction. The heated area is also more detectable.

Conclusion

To optimise the performance of aircraft engines, contact between blades and covers of casing must be analysed, especially at very high speed. Furthermore, numerical simulations require experimental data to validate the modelling. This work presents an experimental device able to attempt at first speeds of interaction of 107 m/s. Orthogonal conditions of interaction are used to measure generated forces. The measurement forces are proportional to the interaction depth and are superimposed on forces due to the friction near the interface of two materials. The increase of the interaction speed imposes a force increase more marked for the high interaction depths. The sharp edge of one interaction material influences essentially the intensity of the friction component of the force. A blunted edge tends to increase the interaction forces. These first results show the influence of the blade tip geometry on the applied interaction load and must be considered to design the blade.

High-speed video recordings allowed to define the nature of shavings produced at high velocities. They are mainly consisted in fine particles. The measurement of the temperature fields during interaction is limited to speeds around 5 m/s. For this speed range the rise of temperature does not seem significant.

Fig. 1 .

 1 Fig. 1. Experimental device with high-speed launcher.

Fig. 3 .

 3 Fig. 3. Detail of the different components used to simulate high-speed interaction.

Fig. 4 .

 4 Fig. 4. Responses of load sensors with strain gauges (a) and KISTLER piezoelectric sensor (b) in terms of quasi-static load.

Fig. 7 .

 7 Fig. 7. Detail of the contact localized under the edge of specimen B.

Fig. 8 .

 8 Fig. 8. Interaction forces as function of translation velocity. Depth of interaction t 1 = 0.05, 0.1, 0.5 mm, edge radius r = 0.001 mm.

Fig. 9 .

 9 Fig. 9. Interaction forces for two edge radii: r = 0.001 and 0.03 mm.

Fig. 10 .

 10 Fig. 10. Photograph of the interaction area in static position. The abradable material (specimen A) is in position near the titanium alloy (specimen B).

Fig. 11 .

 11 Fig. 11. Real time photographs of interaction. (a) Interaction speed V = 46 m/s; depth t 1 = 0.5 mm. (b) Interaction speed V = 56 m/s; depth t 1 = 0.5 mm.

Fig. 12 .

 12 Fig. 12. Temperature field during the interaction abradable/titanium at a speed of 5 m/s.

Table 1

 1 Chemical composition of abradable material M601

		Element (%)
		Si	Fe	C	Impurities	Al	Polyester
	Aluminum-silicon alloy and polyester
	Min.	6.50	-	-	-	51	The rest
	Max.	8.50	0.5	0.2	0.35	55	The rest

Acknowledgment

The authors are grateful for the support provided by SNECMA.


Li Wang

Ali Mohammad-Djafari

Nicolas Gac

X-RAY COMPUTED TOMOGRAPHY SIMULTANEOUS IMAGE RECONSTRUCTION AND CONTOUR DETECTION USING A HIERARCHICAL MARKOVIAN MODEL

Keywords: Computed Tomography, Non Destructive Testing, Simultaneous reconstruction and contour detection, Markovian Model, and Generalized Student-t distribution

The 3D X-ray Computed Tomography (CT) is used in many domains. In medical imaging and industrial Non Destructive Testing (NDT) applications, this technique becomes of great interest. In these applications, very often, we need not only to reconstruct the image, but also to detect the contours between the homogeneous regions of the piecewise continuous image. Generally, contours are obtained by a post processing from the reconstructed image. In this paper, we propose a method to estimate image and contour simultaneously. For this we use the Bayesian approach with a prior model in which the relationship between the image and its contour is considered by using a hierarchical Markovian model, and use a sparsity enforcing prior model for the contours. This proposed method can be used for reconstructions when the image is piecewise continuous. The simulation results are compared with some state of the art methods, and they show the efficiency of simultaneously reconstructing and edge detecting by using proposed method.

INTRODUCTION

X-ray CT is used in many domains. The Radon transform is one of the commonly used model to describe it.

X-ray Computed Tomography

In the imaging systems, objects are often observed by using different techniques: X-rays [START_REF] Kalender | X-ray Computed Tomography[END_REF], ultrasound, microwave, infra-red, etc. Each observation process provides different information on the object under study. In X-ray CT, the detected data provides information on the density distribution of the material inside the object. The CT image plays an important role in medical and industrial Non Destructive Testing (NDT) applications. Its crosssectional images are used for diagnostic and therapeutic purposes in various medical disciplines. Industrial CT scanning has been used in many areas of industry for internal inspection of components. In the 3D CT problem, the Graphics Processing Unit (GPU), presented in details in [START_REF] Gac | Multi GPU parallelization of 3D Bayesian CT algorithm and its application on real foam reonconstruction with incomplete data set[END_REF], is always used to accelerate the computation.

In the work of Mumford and Shah in [START_REF] Mumford | Optimal approximations by piecewise smooth functions and associated variational problems[END_REF], they introduced the method which estimate the image and at the same time detect the segmentations and the contours. They estimate the variables by regularized optimization, in which the l 1 penalty is used to enforce sparsity. The inconvenient in their work is to chose the parameter values for all the regularization terms. In our work, we use a Bayesian method to estimate the image and the contours, estimating all the parameters at the same time. In this Bayesian method, a hierarchical Markovian prior model is defined, which combines the image and the contour informations. We use the heavy tailed distribution to enforce the sparsity of the contour. By using the contour information in the reconstruction of image, the edges of image are preserved. Comparing with the work of Mumford and Shah, our proposed method considers the relationships between (1) contour-projection gradient, (2) contour-image, (3) image-projection and (4) projection-projection gradient. The sparseness of contour is enforced and controlled by using a generalized Student-t distribution.

Forward model: Radon Transform

One of the most commonly used mathematical modelling of CT imaging is the Radon Transform (RT), presented with details in [START_REF] Stanley R Deans | The Radon Transform and some of its applications[END_REF]. It can be written as:

g(r, φ) = L r,φ f (x, y)dl (1)
where f (x, y) represents the image, g the detected projection data, r the perpendicular length from center point of coordinate and φ the considered X ray angle. L r,φ is the length of ray (r, φ) passing through the image. One property that we consider in this work is:

f (x, y) → RT → g(r, φ), f (x, y) → RT → g(r, φ) (2)
where

f (x, y) = ∂ 2 ∂x 2 + ∂ 2 ∂y 2 f (x, y), g(r, φ) = ∂ 2 ∂r 2 g(r, φ) (3)
In the discretized form, the Radon Transform in Equation ( 1) can be written as:

g k = H k f + k or g = Hf + (4)
where g k is a vector representing the projection at angle φ k and g ∈ R M ×1 represents all these vectors putting together.

H k and H ∈ R M ×N are their corresponding linear forward projection operators. Noted by f the discretized f (x, y) and by gk discretized g(r, φ), Eq.( 2) and Eq.( 3) therefore are written as:

g = Hf + , g = H f + ¨ g = D 1 g , f = D 2 f (5)
where D 1 represents the 1-dimensional second order derivation of g which can be implemented as the following convolution operation:

D 1 g = [-1 2 -1] * g,
and D 2 represents the discretized Laplacian operator which can be implemented by the following 2D convolution:

D 2 f =   0 -1 0 -1 4 -1 0 -1 0   * f .
The classical image reconstruction aims at estimating f with g. In this paper, we propose to use the forward relations in [START_REF] Gu | Bayesian sparse estimation using double lomax priors[END_REF] to simultaneously reconstruct the image f and detect the contours f by using a Bayesian approach. For this, we use a hierarchical Markovian prior model to define the image f . Its hidden contour variables q, which is related to f , is used in this prior model. A sparsity enforcing Generalized Student-t prior distribution is used to define the contours.

PROPOSED METHOD OF RECONSTRUCTION OF OBJECT WHILE CONSIDERING CONTOURS

Forward model

Starting from the forward model in ( 5) and assuming to be iid Gaussian distributed with variance v , we quickly get the likelihood:

p(g|f , v ) = N (g|Hf , v I) p(g| f , v ) = N (g|H f , 4v I). (6)
It is not difficult to demonstrate that the variance of ¨ , which is the second order derivation of , equals to 4v . The variance v is supposed to be Inverse-Gamma distributed, so that the value is always positive and approaching to zero.

Sparsity enforcing Student-t prior model for f

When considering the piecewise-continuous image, its Laplacian is sparse with most of the part approaching to zero except the contours. As f is sparse, we use a sparsity enforcing prior to define it. Typically three categories of such distributions are used to enforce sparseness: the Generalized Gaussian (GG) distributions, the Gaussian mixture distributions and the heavy-tailed distributions. The generalized Student-t distribution is heavy-tailed, with the convenient that it can be expressed as the marginal of a Normal-Inverse Gamma distribution:

St(x|α, β) = N (x|0, z)IG(z|α, β) dz. (7)
With this property, we can set the following prior model for f :

p( f |v b ) = N ( f |0, v b I) p(v b |α, β) = IG(v b |α, β). ( 8
)
where α, β are two hyper-parameters of generalized Studentt distribution. With these two parameters, comparing with one-parametered standard Student-t distribution, the sparsity properties can be controlled more precisely.

The image-contour hierarchical Markovian model

By using Bayesian approach, prior informations are modeled as probabilistic distributions. Many prior models which are suitable for a Bayesian piecewise continuous CT image reconstruction method are discussed in the references [START_REF] Gu | Bayesian sparse estimation using double lomax priors[END_REF][START_REF] Ting | Sparse image reconstruction for molecular imaging[END_REF][START_REF] Chen | Nonlocal prior Bayesian tomographic reconstruction[END_REF][START_REF] Molina | Variational posterior distribution approximation in bayesian emission tomography reconstruction using a Gamma mixture prior[END_REF].

In this paper we use a contour-depending Markovian prior model. Let note by f (r) the value of the pixel at the position r and by N r all the four neighbours pixels r . We use the following hierarchical contour-depending Markovian model to define f :

p(f (r)|f -r , q, v a ) = N (f (r)| 1 r ∈Nr (1 -q r ) r ∈Nr (1 -q r )f r , v a ) ∝ v -1 2 a exp    - 1 2 v -1 a f (r) -r ∈Nr (1 -q r )f r r ∈Nr (1 -q r ) 2    (9)
where f -r represents all the elements in vector f except f (r). f r represents the four neighbour pixels of f r and v a is the variance of the normal distribution. q is the normalized f , with the definition: q = f -min f / max f -min f . This distribution describes the fact that: the influence of f r on f r depends on the gradient values at the positions r . For each neighbour pixel, when its gradient value is bigger, its influence on the current pixel is smaller. An example is given in Fig. [START_REF] Kalender | X-ray Computed Tomography[END_REF]. The pixel f r in this figure belongs to a Gaussian prior distribution:

p(f r |f r , q r , v a ) = N f r | 1 2 f r 2 + f r 3 , v a .

The model of the system

With the previous propositions, we obtain a hierarchical model which reconstructs the image and contours simultaneously. f belongs to a Gaussian prior distribution which depends on f , and f belongs to a Generalized Student-t distribution which enforces the sparsity. By combining the proposed prior distributions of f , of f and of all the other hyper parameters, we obtain the statistical distribution model of the system:

                   p(g|f , v ) = N (g|Hf , v I) p(g| f , v ) = N (g|H f , 4v I) p(f |q, v a ) = N (f | 1 r (1-q r ) r (1 -q r )f r , v a I) p( f |v b ) = N ( f |0, Σ b ) where Σ b = diag [v b ] p(v b |α b0 , β b0 ) = j IG(v bj |α b0 , β b0 ) p(v |α 0 , β 0 ) = IG(v |α 0 , β 0 ) (10)
where q is defined below Eq.( 9).

The Bayes method

The Bayesian inference [START_REF] Mark | Basic bayesian methods[END_REF] defines the relationship between prior and posterior distributions by:

p(f, θ|g) = p(g|f, θ)p(f |θ)p(θ) p(g) . (11)
By using the Bayes rule, with all the prior distributions in Equation [START_REF] Qi | Resolution and noise properties of MAP reconstruction for fully 3-D PET[END_REF], the posterior distribution is obtained:

p(f , f , v , v b |g, g) ∝ p(g|f , v )p(g| f , v )• p(f |q)p( f |v b )p(v b |α b0 , β b0 )p(v |α 0 , β 0 ) (12)
Here the conditional probabilities: g knowing f and g knowing f , are independent.

THE ESTIMATION OF VARIABLES

A most commonly used estimation method in Bayesian approach is the Joint Maximum A Posterior (JMAP), presented in [START_REF] Qi | Resolution and noise properties of MAP reconstruction for fully 3-D PET[END_REF], which estimate the variables and parameters as the mode of the posterior distribution, by maximizing the posterior distribution:

f , f , v , v b JM AP = arg max f , f ,v ,v b p f , f , v , v b |g, g (13)

The iteration updating rules

By maximizing the posterior distribution in Eq.( 13), we get the updating equations. For f and f , we use the descent gradient algorithm in order to avoid the calculation of the big size matrix inversion, H t H -1 :

iter : f (k+1) = f (k) -γ (k) 1 ∇J f (k) (14) iter : f (k+1) g = f (k) g -γ (k) 2 ∇J f (k) g (15) f = f g + D2 f /2 (16) v = β 0 + 1 2 g -H f 2 + 1 8 g -H f 2 α 0 + M + 1 (17) v b j = β b 0 + 1 2 f 2 j α b 0 + 3 2 (18)
where

                                                     γ1 = ∇J ( f (k) ) 2 v -1 H∇J ( f (k) ) 2 +4v -1 a ∇J ( f (k)
)

J ( f ) = 1 2 v -1 g -H f 2 + 1 2 v -1 a j fj - 1 r prj r prj frj 2 ∇J ( f ) = -v -1 H t g -H f + v -1 a f -1 p * hr p f * hr γ2 = ∇J ( f ) 2 (4 v ) -1 H∇J ( f ) 2 + Y b ∇J ( f ) 2 where Y b = Σ -1 2 b J ( f ) = 1 2 (4 v ) -1 g -H f 2 + 1 2 Y b f 2 ∇J ( f ) = -(4 v ) -1 H t g -H f + Y t b Y b f (19)
where p = 1-q, p f = p.× f and h r = [0 1 0; 1 0 1; 0 1 0] is a filter which calculates the sum of four neighbour pixels. Fig. [START_REF] Gac | Multi GPU parallelization of 3D Bayesian CT algorithm and its application on real foam reonconstruction with incomplete data set[END_REF] shows the directed acyclic graph (DAG) which presents the principal relationship of the updating of all the variables in the iterative updating rule in Eq.( 14)-Eq.( 19).

g -Eq.14 -f - -D 2 - f f ? + h - g -Eq.15 - f g - q 6 f6 - f f ROCC Fig. 2: DAG of proposed ROCC model.

SIMULATION RESULTS

In the simulations, we use ROCC (Reconstruction of Object Considering Contours) to represent the proposed method. Two other hierarchical models are also used in the simulations: the model with f belonging to a Gaussian distribution (GaussianF), and belonging to a Student-t distribution (Stu-dentF). Fig.

(3) compares the performance of different methods applying on the Shepp Logan image of size 128 × 128. The Least Square(LS) method, Quadratic Regularization(QR) method, presented in [START_REF] Elad | Analysis versus synthesis in signal priors[END_REF][START_REF] Rit | Comparison of analytic and algebraic methods for motioncompensated cone-beam CT reconstruction of the thorax[END_REF], and the two hierarchical modelled methods are compared with the proposed ROCC method. The LS and QR methods performs poorly for this ill-posed inverse problem. The GaussianF method and StudentF method performs well when considering only the quality of the reconstructed image and the projections, but it can be seen that the ROCC method works better on the image reconstruction and at the same time gives more precise reconstructed contours, in which the contour is not used. Logan image with 180 projections. These three figures, from top to bottom, show respectively the relative errors of the image ff 2 / f 2 , the projections gg 2 / g 2 and the contours ff 2 / f 2 its contours. We see from the reconstruction results that the small details of the image are retained, while the contours are obtained simultaneously. The performance of different methods have also be compared by reconstructing the Walnut image of size 256 × 256. The comparison of the relative error of the reconstructed image posed ROCC method always works better than the other state of the art methods.

CONCLUSION

In this paper we presented the CT image reconstruction by using a hierarchical Markovian model. This model preserves the edges of image, and reconstructs simultaneously the image and contours. Simulation results and comparisons with the state of the art methods shows that our proposed ROCC method performs better both on reconstruction of object and contours. This model can also be solved by using the Posterior Mean via Variational Bayesian Approach (VBA) method. Simulations in 3D application will also be done in our future work.

Fig. 1 :

 1 Fig. 1: Example of the contour-depending Markovian structured prior model.

Fig.( 4 )Fig. 3 :

 43 Fig. 3: Comparison of all the methods applying with the Shepp

Fig. 4 :Fig. 5 :

 45 Fig. 4: The reconstruction results for the Shepp Logan image. Three images correspond to (a)the original image, (b)the reconstructed image and (c)the reconstructed contours.

Fig. 6 :

 6 Fig. 6: The reconstructed results for the Walnut image. Three images correspond to (a)the original image, (b)the reconstructed image and (c)the reconstructed contours.


Li Wang
email: li.wang@lss.supelec.fr

Ali Mohammad-Djafari

Nicolas Gac
email: nicolas.gac@lss.supelec.fr

X-ray Computed Tomography using a sparsity enforcing prior model based on Haar transformation in a Bayesian framework

Keywords: Computed Tomography (CT), Bayesian Approach, Hierarchical Model, Generalized Student-t distribution, Joint Maximum A Posterior (JMAP)

published or not. The documents may come

X-ray Computed Tomography using a sparsity enforcing prior model based on Haar transformation in a Bayesian framework

Introduction

Computed Tomography is a technique widely used to achieve medical decisions [START_REF] Abrams | Medical implications of computed tomography (CAT scanning)[END_REF][START_REF] Arridge | Optical tomography in medical imaging[END_REF] and to conduct industrial Non-Destructive Testing (NDT) [START_REF] Hanke | X-ray based methods for non-destructive testing and material characterization[END_REF][START_REF] Halmshaw | Non-destructive testing[END_REF]. In a tomography system, objects (2D or 3D) are observed by using different techniques: X-rays [START_REF] Kalender | X-ray Computed Tomography[END_REF], ultrasound [START_REF] Natterer | A propagation-backpropagation method for ultrasound tomography[END_REF], microwave [START_REF] Bolomey | Microwave tomography: from theory to practical imaging systems[END_REF], infra-red [START_REF] Nelson | Infra-red tomography of port-wine-stain blood vessels in human skin[END_REF], etc. Each observation process provides different type of information on the object under study. In [START_REF] Crowther | The reconstruction of a three-dimensional structure from projections and its application to electron microscopy[END_REF], a pioneering introduction of the reconstruction of a three-dimensional density distribution from its projection data is given. The problem of reconstructing an object from its projections then appeared in various fields of research [START_REF] Kalender | X-ray Computed Tomography[END_REF]. Other reconstruction methods are also proposed, for example the direct reconstruction methods proposed by Bracewell & Riddle [START_REF] Bracewell | Inversion of fan-beam scans in radio astronomy[END_REF] and Berry & Gibbs [START_REF] Berry | The interpretation of optical projections[END_REF].

In X-ray CT, the measured data provides information on the structure of the density inside the object [START_REF] Clackdoyle | Tomographic reconstruction in the 21st century[END_REF]. The main analytical method for CT reconstruction is the Filtered Back-Projection (FBP) [START_REF] Deans | The Radon transform and some of its applications[END_REF][START_REF] Chetih | Tomographic image reconstruction using Filtered Back Projection (FBP) and Algebraic Reconstruction Technique (ART)[END_REF] based on the Radon transform [START_REF] Deans | The Radon transform and some of its applications[END_REF][START_REF] Ludwig | The Radon transform on Euclidean space[END_REF][START_REF] Beylkin | Discrete radon transform[END_REF], with its extension in 3D known as Feldkamp [START_REF] Feldkamp | Practical cone-beam algorithm[END_REF]. It performs fairly well when dealing with problems of reconstruction from sufficient data with a high SNR.

In order to use the least possible radiation dose in medical applications with respect to patients' care, or to spend less time for detection in industrial NDT, reconstruction is expected to be done with ever fewer projections. Then reconstruction becomes an ill-posed inverse problem [START_REF] Oml | Extreme methods for solving ill-posed problems with applications to inverse heat transfer problems[END_REF][START_REF] O'sullivan | A statistical perspective on ill-posed inverse problems[END_REF]. In this case, the conventional analytical reconstruction methods always result artefacts to the reconstructed image. Algebraic methods are used to solve this problem. Here we mention two classical methods: the Algebraic Reconstruction Technique (ART) [START_REF] Gordon | A tutorial on ART (algebraic reconstruction techniques)[END_REF][START_REF] Gordon | Algebraic reconstruction techniques (ART) for three-dimensional electron microscopy and X-ray photography[END_REF] and the Simultaneous Iterative Reconstruction Technique (SIRT) [START_REF] Trampert | Simultaneous iterative reconstruction technique: physical interpretation based on the generalized least squares solution[END_REF][START_REF] Bangliang | The use of simultaneous iterative reconstruction technique for electrical capacitance tomography[END_REF]. These methods as well as many other related algorithms consider the discretized model: g = Hf , where f ∈ R N ×1 represents the image to be reconstructed and g ∈ R M ×1 represents the measured data. Matrix H ∈ R M ×N is the linear projection operator, based either on a line integral or strip or fan beam or cone beam forward operator.

As in general the above-mentioned equation system is underdetermined, there is a need for the regularization methods. These methods approach the final result by considering the discretized model:

g = Hf + , (1)
where ∈ R M ×1 corresponds to the noise and the uncertainties of the system. These methods minimize a criterion consisting of two terms: a discrepancy of the observed data to the output of the forward model, very often g -Hf 2 , and a penalty term Ω(f ). The solution is then defined as:

f = arg min f g -Hf 2 + λΩ(f ) , (2)
where λ is called the regularization parameter. Different regularization methods differ on the choice of Ω(f ). When Ω(f ) = 0, it refers to the Least-Squares (LS) method [START_REF] Ben-Israel | Generalized inverses: theory and applications[END_REF], with the drawback that the obtained solution may be too sensitive to the noise owing to the ill-posedness of the problem and the illconditioning of the operator H. When Ω(f ) = Df 2 , the Quadratic Regularization (QR) method [START_REF] Rit | Comparison of analytic and algebraic methods for motion-compensated cone-beam CT reconstruction of the thorax[END_REF] is implied, where the linear operator D is often a derivation operator. The famous Total Variation (TV) method [START_REF] Chambolle | Image recovery via total variation minimization and related problems[END_REF][START_REF] Goldstein | The Split Bregman method for L1-regularized problems[END_REF] is defined by restricting Ω(f ) = D T V f 1 where D T V is the gradient operator. By using the l 1 norm, the sparsity of the D T V f term is enforced. The appearance of the non-differentiable l 1 term leads to difficulties for the implementation of the optimization algorithm. A large number of methods have been proposed to solve this l 1 norm optimization problem, for example the Newton's method [START_REF] Chan | A nonlinear primal-dual method for Total Variation-based image restoration[END_REF], the Split Bregman method [START_REF] Goldstein | The Split Bregman method for L1-regularized problems[END_REF], etc. Many other efforts have been made to develop fast convergence algorithms to reach the solution of Eq. [START_REF] Arridge | Optical tomography in medical imaging[END_REF]. We may mention here ADMM [START_REF] Wahlberg | An ADMM algorithm for a class of total variation regularized estimation problems[END_REF], ISTA [START_REF] Rolfs | Iterative thresholding algorithm for sparse inverse covariance estimation[END_REF], FISTA [START_REF] Tao | Local Linear Convergence of ISTA and FISTA on the LASSO Problem[END_REF], etc. Other constraints exit to enforce the compact regions, for example the Mumford-Shah functional [START_REF] Mumford | Optimal approximations by piecewise smooth functions and associated variational problems[END_REF].

It should be noted that in the above-mentioned regularization methods, there is a parameter, λ, which controls the trade-off between the data-model adequacy and the regularization penalty. This parameter needs to be fixed. There have been various methods for choosing suitable values for λ, for example the Cross Validation (CV) method and the L-curve method [START_REF] Ramani | Regularization parameter selection for nonlinear iterative image restoration and MRI reconstruction using GCV and SURE-based methods[END_REF][START_REF] Galatsanos | Methods for choosing the regularization parameter and estimating the noise variance in image restoration and their relation[END_REF]. However, the computation for choosing this regularization parameter needs to be done for each simulation and each different set of data, which is expensive and not practical. In some regularization methods, for example the Elastic Net method [START_REF] Zou | Regularization and variable selection via the elastic net[END_REF], there are two or even more such parameters to be fixed.

The data-model adequacy term g -Hf 2 can also take other different expressions, which depend on the noise model. When a Gaussian distribution is used to model the noise, the quadratic expression g -Hf 2 is obtained naturally. In Computed Tomography, in particular in Positrons Emission Tomography when the number of detected photons is very low, the noise may be better modelled by a Poisson distribution. To account for a more precise modelling of the noise, statistical methods have been proposed [START_REF] Wang | Statistical image reconstruction for muon tomography using a Gaussian scale mixture model[END_REF][START_REF] Ramani | Statistical X-ray CT reconstruction using a splitting-based iterative algorithm with orthonormal wavelets[END_REF]. The Maximum Likelihood (ML) methods [START_REF] Redner | Mixture densities, maximum likelihood and the EM algorithm[END_REF] and their different algorithms such as the Expectation Maximization (EM) algorithms [START_REF] Moon | The expectation-maximization algorithm[END_REF][START_REF] Shepp | Maximum likelihood reconstruction for emission tomography[END_REF], the stochastic EM (SEM) [START_REF] Tregouet | A new algorithm for haplotype-based association analysis: the Stochastic-EM algorithm[END_REF], the Ordered subsets-EM (OS-EM) [START_REF] Hudson | Accelerated image reconstruction using ordered subsets of projection data[END_REF] are used commonly in the PET and SPECT reconstruction problems.

In classic ML methods, it is difficult to account for any prior modelling. The Bayesian methods are then widely used [START_REF] Gelman | Bayesian data analysis[END_REF]. The main strategy of the Bayesian inference is to assign a prior probability model to translate our prior knowledge and then to obtain the expression of the posterior distribution. The basic Bayesian formula is:

p(f |g) = p(g|f )p(f ) p(g) , (3)
where p(g|f ) is the likelihood, p(f ) the law a priori and p(f |g) the law a posteriori. Once the posterior distribution is obtained, we can use it for any inference. In practice, however, we may be obliged to use point estimators such as maximum a posteriori (MAP) [START_REF] Qi | Resolution and noise properties of MAP reconstruction for fully 3-D PET[END_REF] or Posterior Mean (PM) [START_REF] Ericson | A note on the posterior mean of a population mean[END_REF]. Using a MAP estimator f = arg max f {p(f |g)} = arg min f {-ln p(f |g)} can make the link between this estimator and almost all the regularization methods. A Gaussian prior for p(f ) results in a Quadratic Regularization method, while a double exponential prior results in a Total Variation regularization method.

When applying the Bayesian approach to real applications, we may need to determine the parameters of the likelihood and the prior models. These parameters, often called hyper-parameters, also have to be estimated in order to achieve unsupervised or semi-supervised methods. Fortunately the Bayesian approach offers this possibility by obtaining the joint posterior probabilistic law:

p(f , θ|g) = p(g|f , θ 1 )p(f |θ 2 )p(θ) p(g) , (4)
where the vector θ = (θ 1 , θ 2 ) contains all the hyper-parameters to be estimated. This joint Posterior law can then be used in different ways to infer on f and θ. One can first marginalize it to obtain p(θ|g) = p(f , θ|g) df and estimate θ = arg max θ {p(θ|g)}, then use it as it was known. Another approach is the joint optimization, which estimates the variable f and hyper-parameters θ simultaneously. The Variational Bayesian Approaching (VBA) method [START_REF] Fox | A tutorial on variational Bayesian inference[END_REF][START_REF] Tzikas | The variational approximation for Bayesian inference[END_REF] is often used during the point estimations.

Often, in order to distinguish the details of an object, a high resolution image is needed. In industrial applications, especially for the NDT of a large-sized object, the size of data becomes critical, and the projection and back-projection operators in CT become expensive. It is then necessary to account for the computational aspects, for example the use of the GPU processor [START_REF] Chen | A hybrid architecture for compressive sensing 3-D CT reconstruction[END_REF][START_REF] Noël | GPU-based cone beam computed tomography[END_REF][START_REF] Gac | High speed 3D tomography on CPU, GPU, and FPGA[END_REF].

In the previous work [START_REF] Wang | Computed tomography reconstruction based on a hierarchical model and variational Bayesian method[END_REF], we proposed using the Bayesian method via a synthesis model f = Dz, in which the Haar transformation coefficient z of the image is first estimated, and then the final image reconstruction result is obtained from the inversion of the Haar transformation of the coefficients by f = D z. When using a Double Exponential (DE) prior model for z and the MAP estimator, the problem becomes equivalent to the optimization of J (z) = g -HDz 2 + λ z 1 which has been used very often in l 1 regularization methods. The specificity of our work was to use a generalized Student-t model in place of DE. In this paper, we present a method which estimate the image and the transformation coefficient simultaneously. The generalization of Student-t distribution is used to enforce the sparsity of the Haar transformation coefficient of the image. In this way, the edges of the reconstructed image is preserved compared with the previous work. From the reconstructed Haar transformation coefficient, we also obtain the high-frequency component (the edge) information of the image.

The rest of this paper is organized as follows. First, we present the Multilevel Haar transformation which is used in the proposed model. Then some sparsity-enforcing distributions and a generalization of Student-t distribution is presented. The proposed Bayesian method is then introduced with details. Finally, we give our simulation results and draw our conclusions.

Notation: Column vectors are represented as boldface lower-case letters, e.g., x, and matrices as boldface upper-case letters, e.g., X, while the i-th component of vector x is denoted by x i the ij-th element of matrix X by X ij . Moreover, (•) T denotes transposition, • 1 stands for the l 1 -norm, • stands for the standard l 2 -norm, |•| denotes the determinant of a matrix or absolute value in the case of a scalar, N (•) is the Gaussian distribution, IG(•) is the Inverse-Gamma (IG) distribution, diag [x] denotes a diagonal matrix whose diagonal entries are the elements of x.

The Hierarchical Haar transform-Based Method (HHBM)

In this paper, we propose a HHBM method, which considers to reconstruct X-ray CT images by using a Bayesian method. This method can be used to solve the CT reconstruction problems for piecewise continuous images consisting of several homogeneous blocs. By using this proposed method, we impose a hidden variable: the multilevel Haar transformation coefficient of the image. This transformation coefficient has a sparse property, which is defined in this method by a generalized Student-t distribution. With the simultaneous estimations of the coefficient and the image, the edges of the image are preserved.

Multilevel Haar Transformation

In a typical X-ray CT problem, the object under consideration consists of a finite number of homogeneous materials, leading to the fact that the images to be reconstructed are piecewise continuous. In each homogeneous bloc which corresponds to one material, the gradient of the image is zero. There are non-zeros values only at the edges between different blocs. In order to take into account this property, we use a sparsity-enforcing law to define a sparse transformation coefficient of the object.

The piece-wise continuous objects [START_REF] Elad | Image denoising via sparse and redundant representations over learned dictionaries[END_REF][START_REF] Chen | Hyperspectral image classification using dictionary-based sparse representation[END_REF] can be represented by different kinds of sparse transformation coefficients. In this paper, we use the multilevel Haar transformation, shown in Fig. [START_REF] Abrams | Medical implications of computed tomography (CAT scanning)[END_REF], with the following properties:

• It is an orthogonal transformation [START_REF] Flury | An algorithm for simultaneous orthogonal transformation of several positive definite symmetric matrices to nearly diagonal form[END_REF][START_REF] Montgomery | Topological transformation groups[END_REF], meaning that if we note by f = Dz we have z =

D -1 f = D T f , with the property D T = D -1 .
• It has a multilevel structure that, after the transformation, only the low-pass component, which corresponds to the top-left small part in the coefficient image shown in Fig. [START_REF] Abrams | Medical implications of computed tomography (CAT scanning)[END_REF], is not sparse.

• Both the forward and inverse operators of the transformation can be implemented easily for 2D and 3D cases. For an object of size N , the cost of computation is only 2 × o(log N ) or 3 × o(log N ) for 2D or 3D cases respectively.

The transformation operator can be represented by:

f = Dz + ξ, (5)
where ξ represents the noise of this transformation.

For the multilevel Haar transformation, we define the transform level by l: f

(l) H = H l (f ).
The definition of this operator is:

f (i) H = H i (f ) = Haar f (i-1) H 1 : N 2 i-1 , 1 : N 2 i-1 , (6) f (1)
H = H 1 (f ) = Haar (f (1 : N, 1 : N )) , (7)
where Haar(•) represents the Haar transformation [START_REF] Ahmed | Discrete cosine transform[END_REF].

Generalized Student-t Distribution

The sparsity of the multilevel Haar transformation coefficient z need to be defined by a statistical model. Typically the sparse property of a statistical variable is enforced by using three kinds of distributions: the Generalized Gaussian distributions [START_REF] Nadarajah | A generalized normal distribution[END_REF], the Gaussian Mixture distributions [START_REF] Rasmussen | The infinite Gaussian mixture model[END_REF][START_REF] Aristophanous | A Gaussian mixture model for definition of lung tumor volumes in positron emission tomography[END_REF] and the heavy tailed distributions [START_REF] Klebanov | Heavy tailed distributions[END_REF][START_REF] Gomes | Heavy-tailed phenomena in satisfiability and constraint satisfaction problems[END_REF]. The Student-t distribution [START_REF] Lange | Robust statistical modeling using the t distribution[END_REF][START_REF] Tzikas | Variational Bayesian sparse kernel-based blind image deconvolution with Student's-t priors[END_REF] is heavy tailed. We present here a generalization of this distribution [START_REF] Dumitru | A Bayesian approach for periodic components estimation for chronobiological signals[END_REF], which has some good properties for the definition of the prior model. The standard Student-t distribution has the probability density function (pdf):

St(x|v) = Γ v+1 2 √ vπΓ v 2 1 + x 2 v -v+1 2 , (8)
of which the variance depends only on the degree of freedom v:

Var [x] =      v v-2 for v > 2 ∞ for 1 < v < 2 undefined otherwise . ( 9
)
From the definition we easily determine that there is a limit for the variance of the Student-t distribution:

Var [x] = v v-2 = 1 + 2 v-2 > 1.
This limit leads to the consequence that this heavy-tailed distribution cannot have a small variance, therefore the sparsity could not be intensively enforced.

We introduce the generalized Student-t distribution (Stg). It comes from the Normal-Gamma marginalization property:

St g (x|α, β) = N (x|0, 1/v) G (v|α, β) dv or the Normal-Inverse Gamma marginal- ization property: St g (x|α, β) = N (x|0, v) IG (v|α, β) dv of the Student-t distribution
Theorem 2.1. For two random variables x and v following a bi-variate Normal-Inverse Gamma distribution: (x, v) ∼ N (x|0, v) IG (v|α, β), the marginal distribution corresponding to x is a Student-t like distribution, generalizing the classical form, denoted St g , depending on the two parameters α and β and having the following pdf:

St g (x|α, β) = N (x|0, v) IG (v|α, β) dv = Γ α + 1 2 √ 2πβΓ (α) 1 + x 2 2β -(α+ 1
2 ) .

(

) 10
This generalization of the Student-t distribution adds a supplementary parameter to the standard one. More possibilities of the distribution shape are therefore available. The shape of the peak and tail of the distribution can be controlled by changing the values of these two parameters. We can either get a super heavy-tailed distribution with a super narrow peak, or a distribution that approaches the Normal distribution, by setting different pair of parameters α and β.

Remark 2.2. If x belong to a generalized Student-t distribution, then it has the expectation and variance:

E [x] = 0 for α > 1 2 undef ined otherwise , (11)
and

Var [x] =      β α-1 for α > 1 ∞ for 1 2 < α ≤ 1 undef ined for 0 < α ≤ 1 2 . ( 12
)
With this property, we can use the St g distribution to define the hidden variables, the multilevel Haar transformation coefficient in the HHBM method, or other sparse representations of the image.

In this paper, we use this generalization of Student-t distribution to define the sparse coefficient. First of all, by commanding the two hyper-parameters, we can control the sparsity rate of the coefficient. What's more, thanks to the fact that it is a marginalization of two conjugate distributions, it is more convenient for us to apply different estimation methods, for example the Posterior Mean (PM) via Variational Bayesian Approximation (VBA).

Hierarchical Bayesian Model of HHBM

The Bayesian method attempts to solve the reconstruction problem statistically based on the forward model presented in Eq. [START_REF] Abrams | Medical implications of computed tomography (CAT scanning)[END_REF]. By using the Bayesian method, first of all we define the probabilistic model for all the prior information.

The additive noise in this model is defined as belonging to a Gaussian distribution, with mean equal to zero and variance v :

p( |v ) = N ( |0, V ) , (13)
where V = diag [v ] is a diagonal matrix. According to the linear forward model, the distribution of the additive noise provides the likelihood:

p(g|f , v ) = N (g|Hf , V ). ( 14
)
For the multilevel Haar transformation of the image, the likelihood corresponding to the forward model in Eq.( 5) is given below:

p(f |z, v ξ ) = N (f |Dz, V ξ ), (15)
where D represents the inverse Multilevel Haar transformation operator, and the additive noise ξ is considered to belong to a Gaussian distribution, with mean equal to zero and variance v ξ . The definition of the covariance matrix is:

V ξ = diag [v ξ ]. The vector z = [z 1 , z 2 , • • • , z N ],
as mentioned above, is sparse. As discussed in Section (2.2), the generalized Student-t distribution is used to define z in order to enforce the sparsity property. By using the Normal-Inverse Gamma property shown in Eq.( 10), the prior distribution for z can be realized by a hierarchical structured Normal-Inverse Gamma distribution, where z belongs to a Normal distribution, the variance of which belongs to an Inverse Gamma distribution depending on two hyper-parameters α z 0 and β z 0 :

p(z|v z ) = N (z|0, V z ), (16)
p(v z |α z 0 , β z 0 ) = N j=1 IG(v z j |α z 0 , β z 0 ), (17)
where the elements of vector v z are supposed to be i.i.d. Moreover, we need to decide the probabilistic model for the variances of the two noises: v and v ξ in this unsupervised model. Knowing that the variance is always positive, and that typically the value of the variance is small, we choose the Inverse Gamma distribution to model them:

p(v |α 0 , β 0 ) = M i=1 IG(v i |α 0 , β 0 ), (18)
p(v ξ |α ξ 0 , β ξ 0 ) = N j=1 IG(v ξ j |α ξ 0 , β ξ 0 ), (19)
where the elements in each vector, v and v ξ , are supposed to be i.i.d. α 0 , β 0 , α ξ 0 and β ξ 0 are hyperparameters for these distributions respectively.

As long as we have defined the models for all the prior distributions, the model of the hierarchical system containing all the variables, parameters and hyper-parameters is:

p(g|f , v ) = N (g|Hf , V ) ∝ |V | -1 2 exp - 1 2 (g -Hf ) T V -1 (g -Hf ) (20)
p(f |z, v ξ ) = N (f |Dz, V ξ ) ∝ |V ξ | -1 2 exp - 1 2 (f -Dz) T V -1 ξ (f -Dz) (21)
p(z|v z ) = N (z|0, V z ) ∝ |V z | -1 2 exp - 1 2 z T V -1 z z (22) p(v z |α z0 , β z0 ) = N j=1 IG(v zj |α z0 , β z0 ) ∝ N j=1 v -(αz 0 +1) zj exp -β z0 v -1 zj (23) p(v |α 0 , β 0 ) = M i=1 IG(v i |α 0 , β 0 ) ∝ M i=1 v -(α 0 +1) i exp -β 0 v -1 i (24) p(v ξ |α ξ0 , β ξ0 ) = N j=1 IG(v ξj |α ξ0 , β ξ0 ) ∝ N j=1 v -(α ξ 0 +1) ξj exp -β ξ0 v -1 ξj (25)
The corresponding directed acyclic graph (DAG) of the proposed model is shown in Fig. [START_REF] Arridge | Optical tomography in medical imaging[END_REF]. The posterior distribution is then obtained from Equations Eq.( 20)-Eq.( 25) via unsupervised Bayesian formula:

α z 0 , β z 0 ? v z ? z ? D α 0 , β 0 ? v ? α ξ 0 , β ξ 0 ? v ξ ? ξ -f ? H g
p(f , z, v , v ξ , v z |g) ∝ p(g|f , v )p(f |z, v ξ )p(z|v z )p(v )p(v ξ )p(v z ) = |V | -1 2 exp - 1 2 (g -Hf ) T V -1 (g -Hf ) |V ξ | -1 2 exp - 1 2 (f -Dz) T V -1 ξ (f -Dz) |V z | -1 2 exp - 1 2 z T V -1 z z N j=1 v -(αz 0 +1) z j exp -β z 0 v -1 z j M i=1 v -(α 0 +1) i exp -β 0 v -1 i N j=1 v -(α ξ 0 +1) ξ j exp -β ξ 0 v -1 ξ j (26)
Once the posterior distribution has been obtained, estimation methods, for example the MAP algorithm [START_REF] Qi | Resolution and noise properties of MAP reconstruction for fully 3-D PET[END_REF] or the posterior mean (PM) via variational Bayesian approach (VBA) [START_REF] Fox | A tutorial on variational Bayesian inference[END_REF][START_REF] Tzikas | The variational approximation for Bayesian inference[END_REF] method, are then used to approach the reconstruction result. In this paper we use the Joint MAP method which iteratively estimates all the variables and parameters.

Joint Maximum A Posterior (JMAP) Algorithm for HHBM

The JMAP algorithm aims at iteratively and alternately estimating all the variables and parameters by maximizing the posterior distribution, or minimizing the negative logarithm of the posterior law:

( f , z, v , v ξ , v z ) = arg max f ,z,v ,v ξ ,vz {p(f , z, v , v ξ , v z |g)} = arg min f ,z,v ,v ξ ,vz {-ln p(f , z, v , v ξ , v z |g)} . ( 27
)
The criterion of the minimization,

J (f , z, v , v ξ , v z ) = -ln p(f , z, v , v ξ , v z |g)
, is then obtained from Eq.( 26):

J (f , z, v , v ξ , v z ) = 1 2 M i=1 ln v i + 1 2 (g -Hf ) T V -1 (g -Hf ) + 1 2 N j=1 ln v ξ j + 1 2 (f -Dz) T V -1 ξ (f -Dz) + 1 2 N j=1 ln v z j + 1 2 z T V -1 z z + (α z 0 + 1) N j=1 ln v z j + N j=1 β z 0 v -1 z j + (α 0 + 1) M i=1 ln v i + M i=1 β 0 v -1 i + (α ξ 0 + 1) N j=1 ln v ξ j + N j=1 β ξ 0 v -1 ξ j . (28)
By using the JMAP algorithm, we estimate all the variables by minimizing the corresponding terms in the constraint equation Eq.( 28). While considering one variable, we keep all the other variables fixed, equal to the estimated value of the previous iteration.

f (k) = arg min f J (f , z (k-1) , v (k-1) z , v (k-1) , v (k-1) ξ ) , (29)
z (k) = arg min z J ( f (k) , z, v (k-1) z , v (k-1) , v (k-1) ξ ) , (30)
v (k) z = arg min vz J ( f (k) , z (k) , v z , v (k-1) , v (k-1) ξ ) , (31)
v (k) = arg min v J ( f (k) , z (k) , v (k) z , v , v (k-1) ξ ) , (32) v (k)
ξ = arg min v ξ J ( f (k) , z (k) , v (k) z , v (k) , v ξ ) , (33)
where • (k) represents the estimated value of the kth iteration. The variables with the hat symbol mean that in the current step they are considered as knowns. We then get the expressions for the minimization constraints as:

f (k+1) = arg min f 1 2 (g -Hf ) T V (k) -1 (g -Hf ) + 1 2 f -D z (k) T V ξ (k) -1 f -D z (k) , (34)
z (k+1) = arg min z 1 2 f (k+1) -Dz T V ξ (k) -1 f (k+1) -Dz + 1 2 z T V z (k) -1 z , (35)
v (k+1) z = arg min vz    1 2 N j=1 ln vz j + 1 2 z (k+1) T V -1 z z (k+1) + (αz 0 + 1) N j=1 ln vz j + N j=1 βz 0 v -1 z j    , (36)
v (k+1) = arg min v 1 2 M i=1 ln v i + 1 2 g -H f (k+1) T V -1 g -H f (k+1) + (α 0 + 1) M i=1 ln v i + M i=1 β 0 v -1 i , (37)
v (k+1) ξ = arg min v ξ    1 2 N j=1 ln v ξ j + 1 2 f (k+1) -D z (k+1) T V -1 ξ f (k+1) -D z (k+1) + (α ξ 0 + 1) N j=1 ln v ξ j + N j=1 β ξ 0 v -1 ξ j    . (38)

Estimate f and z

For the estimations of the object f and the transformation coefficient z, the updating equations obtained from the minimization in Eq. [START_REF] Galatsanos | Methods for choosing the regularization parameter and estimating the noise variance in image restoration and their relation[END_REF] and Eq.( 35) are:

f = H T V -1 H + V -1 ξ -1 H T V -1 g + V -1 ξ Dz , (39)
z = D T V -1 ξ D + V -1 z -1 D T V -1 ξ f . ( 40
)
As we have already said, in X-ray CT, the size of the object f and the size of its transformation coefficient z are very large. Hence the computation of the matrix inversion in Eq.( 39) and Eq.( 40) is too expensive and even impossible to compute. In this paper we present the gradient descent algorithm for these estimations.

for k = 1 → Iter2 : f (k+1) = f (k) -γ (k) f ∇J ( f (k) ), (41)
for k = 1 → Iter2 : z (k+1) = z (k) -γ (k) z ∇J ( z (k) ), ( 42
)
where Iter2 is the number of iterations for the descent gradient step, ∇J (•) is the corresponding gradient of the constraint equation, and γ (•) is the corresponding descent step length:

∇J (f ) = -H T V -1 (g -Hf ) + V -1 ξ (f -Dz) , (43)
∇J (z) = -D T V -1 ξ (f -Dz) + V -1 z z, (44) γ (k)
f = ∇J ( f (k) ) 2 / Y H∇J ( f (k) ) 2 + Y ξ ∇J ( f (k) ) 2 , (45)
γ (k) z = ∇J ( z (k) ) 2 / Y ξ D∇J ( z (k) ) 2 + Y z ∇J ( z (k) ) 2 . ( 46
)
The descent step length γ f and γ z are obtained by using an optimized step length strategy, see [START_REF] Boyd | Convex optimization[END_REF]. The calculation details are given in Appendix A. The computation algorithms for updating f and z are given in Algo. Calculate ∇J ( f

(k-1)
) according to Eq.(

Update γ Update Calculate ∇J ( z (k-1) ) according to Eq.( 44) Update

f (k) = f (k-1) -γ (k) f ∇J ( f (k-1)
z (k) = z (k-1) -γ (k)
z ∇J ( z (k-1) ) 9: end for

Estimation of the other variables

The iterative updating rules for the other variables are obtained by calculating the minimizations in Eq.( 36)-Eq.( 38) and are given as the following formulas:

v z j = β z 0 + 1 2 z 2 j / (α z 0 + 3/2) , ∀j ∈ [1, N ], (47)
v i = β 0 + 1 2 g i -H f i 2 / (α 0 + 3/2) , ∀i ∈ [1, M ], (48)
v ξ j = β ξ 0 + 1 2 f j -[D z] j 2 / (α ξ 0 + 3/2) , ∀j ∈ [1, N ]. (49)

Summary of the algorithm

Now we have obtained all the updating equations for all the variables. In the computation, we still need to initialize f and z. For the initialization of f , we use the least square estimation:

f ini = arg min f g -Hf 2 (50)
with a very limited number of iterations. Then we initialize z as the l-level Haar transformation coefficient of f ini : z ini = D -1 f ini . The initialization of the other unknowns is obtained from Eq.( 36)-Eq.( 38) by using f ini and z ini . The Algo. [START_REF] Hanke | X-ray based methods for non-destructive testing and material characterization[END_REF] gives the details of the JMAP algorithm for the HHBM method.

Simulations

In the simulations, the 2D Shepp Logan phantom of size 64 × 64 and 256 × 256 are used as the original images. Fig. (3)(a) shows the image of size 256 × 256. It consists of several homogeneous blocs, each of which corresponds to a different material. The projections are obtained from angles uniformly distributed from 0 • to 180 • , no matter how many number of projections we apply.

The comparison is between a) the proposed HHBM method, b) the Quadratic Regularization (QR) method [START_REF] Rit | Comparison of analytic and algebraic methods for motion-compensated cone-beam CT reconstruction of the thorax[END_REF], which is a basic regularization method and c) the Total Variation (TV) method [START_REF] Chambolle | Image recovery via total variation minimization and related problems[END_REF], which performs fairly well when treating piecewise continuous image reconstruction problems. The reconstruction by using the TV method is realized by using the previously mentioned Split Bregman algorithm [START_REF] Goldstein | The Split Bregman method for L1-regularized problems[END_REF]. The gradient descend optimization algorithm is used in the HHBM method as mentioned in Eq.( 41) and Eq.( 42).

Algorithm 3 Summary of the JMAP algorithm for HHBM method 1: Fix parameters α z 0 , β z 0 , α 0 , β 0 , α ξ 0 , β ξ 0 , l 2: Input: H, D, g 3: Output: f , z, v z , v , v ξ 4: Initialization:

5: f (0)
← arg min f g -Hf 2 with a few iterations 6: z (0) ← l-level Haar transformation of f

Measurement of performance

A large number of measurements can be used for the evaluation of the reconstruction performance of the different methods. In this paper, we consider three evaluations: the relative error δ f , the improvement of the signal-to-noise ratio (ISNR), and the peak signal-to-noise ratio (PSNR). They show respectively the quality of the entire reconstructed image, the improvement of the estimation iterations from the initialized image, and the signal-to-noise ratio (SNR) of the result corresponding to the biggest value of the image. The definitions of these evaluations are:

• δ f = f -f 2 f 2
, where f represents the original image and f represents the reconstructed image.

• ISN R = 10 log 10 f -f (0) 2 f -f 2 ,
where f (0) is the initialized image for the iterative algorithms.

• P SN R = 10 log 10 max (f )

2 1 N j (fj-f j )
2 , where N is the length of the vector f .

Simulation results

In the Table ( From the results in these two tables, we find that for the reconstruction of the 64×64 Shepp Logan image, the proposed HHBM method performs quite well compared to the other two state-of-the art methods, in all respects of the evaluation (relative error, ISNR, PSNR). In the reconstruction of the 256 × 256 Shepp Logan image, the proposed HHBM method has a better reconstruction result within 50 iterations when the data is noisy with a high SNR, with a sufficient or insufficient number of projections. When the data is noisy with a low SNR=20dB, the proposed HHBM performs better with respect to the PSNR evaluation, but the TV method surpasses the HHBM method with respect to the relative error value. Discussion of the results will be given after, in the subsection 4.3.

Table (3) shows different metrics of the 256 × 256 Shepp Logan image reconstructed from 128, 64 and 32 projections data set when the following convergence criterion is reached:

δ (k+1) f -δ (k) f < ε, (51)
where ε is a small fixed value, δ

f represents the relative error after k iterations, and in this simulation we set ε = 0.0001. We also estimated the structural similarity (SSIM) metric [START_REF] Wang | Image quality assessment: from error visibility to structural similarity[END_REF] in the comparison, with the definition:

SSIM (x, y) = (2µ x µ y + C 1 ) (2σ xy + C 2 ) µ 2 x + µ 2 y + C 1 σ 2 x + σ 2 y + C 2 , (52)
where C 1 and C 2 are fixed constant, µ is the mean and σ is the standard deviation. We can conclude from the table that the difference of the SSIM of final results by using different methods are not evident. In order to show clearly the reconstructed Shepp Logan phantom, Fig.( 5)- [START_REF] Nelson | Infra-red tomography of port-wine-stain blood vessels in human skin[END_REF] show two zones of the image. The two zones are correspondingly the red and blue areas shown in Fig. (4). In the results, an example of profile is displayed in the images in order to show the change at the contour positions.

The comparison of the reconstructed 2D 256 × 256 images using different methods, with 128 projections and a high SNR=40dB, is shown in Fig. [START_REF] Kalender | X-ray Computed Tomography[END_REF]. In the images, we see that in the reconstructed image using the QR method, the shape of different blocs is not clear. In the result of the TV method, we find that the homogeneous areas of the reconstructed image are very smooth, but compared with the result of the proposed HHBM method, the shape of the blocs is not clear enough. In the HHBM reconstructed image, there are a few artefacts in the homogeneous blocs, but the shape of the image is clearer than these two state-of-the-art methods. In the figures, the blue line is the position of the profile which is displayed in red. The comparison of the reconstructed 2D 256 × 256 images using different methods, with 64 projections and a high SNR=40dB, is shown in Fig. [START_REF] Bolomey | Microwave tomography: from theory to practical imaging systems[END_REF]. In this case, the reconstruction is done by using only a quarter of the projections data. In this case, the TV method begins to have artefacts in the homogeneous areas, and the contours are not very clear. The result of the HHBM method has clearer contours.

The comparison of the reconstructed 2D 256 × 256 images using different methods, with 32 projections and a high SNR=40dB, is shown in Fig. [START_REF] Nelson | Infra-red tomography of port-wine-stain blood vessels in human skin[END_REF]. In this case, the reconstruction is done by using only 1/8 of the projection data. The result of the HHBM method still has clearer contours.

In Fig. [START_REF] Crowther | The reconstruction of a three-dimensional structure from projections and its application to electron microscopy[END_REF], Fig. [START_REF] Bracewell | Inversion of fan-beam scans in radio astronomy[END_REF] and Fig. [START_REF] Berry | The interpretation of optical projections[END_REF], the comparisons of the convergence of the relative error are shown, for the simulations with 128 projections and a high SNR=40dB, with 128 projections and a low SNR=20dB, with 64 projections and a high SNR=40dB, respectively. From these comparisons, we see that by using the HHBM method, the convergence is faster than the other state-of-the-art methods. Within just a few iterations, the reconstructed image already has a very good quality. But in the simulation with a low SNR, the relative error of the final results of HHBM is not very small.

In Fig. [START_REF] Clackdoyle | Tomographic reconstruction in the 21st century[END_REF], we compare the convergence performance of the reconstruction with only 32 projections for a reconstruction of the Shepp Logan phantom of size 256 × 256. In this comparison, we see that the HHBM method performs better than the FBP, QR and TV methods.

The comparison of the profiles of the reconstructed images is given in Fig. see more clearly the details. In these profile figures, we see that by using the HHBM method, the profile is closer to the original profile at the contour positions, thus at the areas where there are major changes. Then we analyze the performance of different methods when using a different number of projections and different SNR. The comparisons are given in Fig. [START_REF] Ludwig | The Radon transform on Euclidean space[END_REF] and Fig. [START_REF] Beylkin | Discrete radon transform[END_REF]. In these two figures, the left one analyzes the influence of the number of projections, and the right one analyzes the influence of the value of SNR. In all these figures, we use the relative error to evaluate the results. The evaluation is measured after 50 iterations for each method and each simulation. We can see from the results that all the methods have better results when dealing with (1) larger numbers of projections, and (2) data with a higher SNR. From the results, we notice that, when reconstructing from the data with a high SNR, the HHBM method is robust, even when the projection number is very small. When the size of phantom is smaller, the TV method performs better when there is a very low SNR (10dB). When the size of phantom becomes bigger, the noise has a significant influence on the TV method.

Discussion

The main drawback of the proposed HHBM method is that the noise in the homogeneous areas cannot be removed efficiently when the data has a low SNR. The problem comes from the hierarchical structured multilevel Haar transformation. In this transformation, the low-frequency component of the image, cor- 17). The problem is that the information of the low-frequency component is more important than the high-frequency component so that in future work, we are considering applying a multi-resolution estimation for the coefficient z. The multi-resolution estimation sets out to estimate at the first step the low-frequency component of z, and then step by step, to estimate the low-frequency component and higher-frequency component.

In the regularization methods, there is always a regularization parameter λ to be fixed. In order to choose a suitable value for this parameter, we compared the relative error of the reconstructed image by using QR and TV methods with λ = [0.1, 0.2, 0.5, [START_REF] Abrams | Medical implications of computed tomography (CAT scanning)[END_REF][START_REF] Bracewell | Inversion of fan-beam scans in radio astronomy[END_REF][START_REF] Gordon | A tutorial on ART (algebraic reconstruction techniques)[END_REF][START_REF] Rolfs | Iterative thresholding algorithm for sparse inverse covariance estimation[END_REF][START_REF] Shepp | Maximum likelihood reconstruction for emission tomography[END_REF][START_REF] Gac | High speed 3D tomography on CPU, GPU, and FPGA[END_REF]100,200] respectively. The results are shown in Fig. [START_REF] Feldkamp | Practical cone-beam algorithm[END_REF].

We can see that the choice of the value of λ plays an important role in the regularization methods. For example in the low SNR cases, it is better to choose a bigger value for λ than the high SNR cases. The cost of choosing such a value is very expensive. In our simulation we tested for only 10 different values and this step takes already much more time than the proposed Bayesian method, needless to say in the real applications where we know nothing about the noise level and the evaluation of the convergence is not appropriate. Hence, the proposed Bayesian method is more appropriate thanks to the fact that such parameters are all estimated automatically during the simulations.

Conclusion and Perspectives

The reconstruction of CT image plays an important role in medical and industrial applications. It is becoming important to reconstruct with an increasingly smaller number of projections. In this paper we proposed a Bayesian method with a hierarchical structured prior model. In this model the multilevel Haar transformation is used as a sparse transformation. A generalized Student-t distribution is used to model the transformation coefficients. The simulation results of the Shepp Logan phantom show good performance of the proposed method compared with the conventional quadratic regularization and Total Variation methods. The proposed method converges much faster than these state-of-the-art methods. In future work, we will look into different sparse transformations, for example the Dual-Tree CWT transformation [START_REF] Selesnick | The dual-tree complex wavelet transform[END_REF] which has a better performance than the Haar transformation. We will also consider applying a 'multi-resolution' estimation strategy for the transformation coefficients. The simulation with 3D data is under execution, with the use of the GPU processor and the ASTRA toolbox [START_REF] Palenstijn | The ASTRA tomography toolbox[END_REF].

A. Expressions for updating step length γ f and γ z A.1. Calculation of parameter γ (k) f

The gradient descent step length for updating f is calculated by minimizing the criterion J f (f -γ f dJ f ):

γ f = arg min γ f J f (f -γ f dJ f ),
where J f (f -γ f dJ f ) is simplified as:

J f (f -γ f dJ f ) = 1 2 g T -f T H T + γ f dJ T f H T V -1 (g -Hf + γ f HdJ f ) + 1 2 f T -γ f dJ T f -z T D T V -1 ξ (f -γ f dJ f -Dz) = 1 2 γ f g T V -1 HdJ f -γ f f T H T V -1 HdJ f + γ f dJ T f H T V -1 g -γ f dJ T f H T V -1 Hf + γ 2 f dJ T f H T V -1 HdJ f + 1 2 -γ f f T V -1 ξ dJ f -γ f dJ T f V -1 ξ f + γ 2 f dJ f V -1 ξ dJ f + γ f dJ T f V -1 ξ Dz +γ f z T D T V -1 ξ dJ f + Cst = γ f dJ T f H T V -1 g -H T V -1 Hf -V -1 ξ f + V -1 ξ Dz + 1 2 γ 2 f Y HdJ f 2 + Y ξ dJ f 2 + Cst.
In Eq.( 43) we have:

∇J (f ) = -H T V -1 (g -Hf ) + V -1 ξ (f -Dz) = -H T V -1 g + H T V -1 Hf + V -1 ξ f -V -1 ξ Dz. So that: J f (f -γ f dJ f ) = arg min γ f -γ f dJ f 2 + 1 2 γ 2 f Y HdJ f 2 + Y ξ dJ f 2 + Cst.
The gradient of the criterion is:

∂J f (f -γ f dJ f ) ∂γ f = -dJ f 2 + γ f Y HdJ f 2 + Y ξ dJ f 2 .
By setting

∂J f (f -γ f dJ f ) ∂γ f | γ f = 0,
we get:

γ f = dJ f 2 Y HdJ f 2 + Y ξ dJ f

Figure 1 :

 1 Figure 1: Multi-level Haar Transformation of the Shepp Logan image. Left: 2-level Haar transformation; right: 3-level Haar transformation.

Figure 2 :

 2 Figure 2: DAG of the proposed HHBM model.

( 1 )Algorithm 1 descent algorithm for f 1 : 2 :: f 3 :) 4 : 5 :

 1112345 and Algo.(2) respectively. Gradient Input: H, D, g, z OutputInitialization: f (0for it = 1 to Iter2 do Set k = it 6:

f

  according to Eq.[START_REF] Ericson | A note on the posterior mean of a population mean[END_REF] 8:

) 9 : 1 : 5 :

 915 end for Algorithm 2 Gradient descent algorithm for z Input: H, D, g, f 2: Output: z 3: Initialization: z (0) 4: for it = 1 to Iter2 do Set k = it 6:

7 :

 7 for it = 1 to Itmax do

  Compute v z , v , v ξ according to Eq.(47)-Eq.[START_REF] Noël | GPU-based cone beam computed tomography[END_REF] 11: end for In the simulations, we have tested with 4 and 5 and 6 levels of transformation, and we found that this transformation level doesn't have big influence on the cost of computation and the final performance. In the simulations, we use a 5-level Haar transformation as the sparse coefficient of the image, shown in Fig.(3)(b). The sinogram of the projection of the 2D 256×256 Shepp Logan phantom with 128 projections and a high SNR=40dB is shown in Fig.(3)(c), each vertical line of which corresponds to the detected data of one projection direction. In Fig.(3)(d), the sinogram of projection with a low SNR=20dB is displayed.

Figure 3 :

 3 Figure 3: (a) original 2D 256 × 256 Shepp Logan phantom, (b) its 5-level Haar transformation, (c) its sinogram of projection with 128 projections and a high SNR=40dB and (d) its sinogram of projection with 128 projection and a low SNR=20dB.

Figure 4 :

 4 Figure 4: The original 256 × 256 Shepp Logan phantom and two zones used to demonstrate the results.

Figure 5 :Figure 6 :

 56 Figure 5: Comparison of zones of reconstructed 256×256 Shepp Logan images by using different methods with 128 projections and a high SNR=40dB by using: (a) QR, (b) TV and (c) proposed HHBM method. The red curve shows the profile at the blue line position.

Figure 7 :Figure 8 :Figure 9 :

 789 Figure 7: Comparison of reconstructed 256 × 256 Shepp Logan images by using different methods with 64 projections and a high SNR=40dB by using: (a) QR, (b) TV and (c) proposed HHBM method. The red curve shows the profile at the blue line position.

Figure 10 :Figure 11 :Figure 12 :

 101112 Figure 10: The performance of convergency of relative error δ f for phantom of size 256 × 256 by using different methods with 128 projections and a low SNR=20dB.

Figure 13 :

 13 Figure 13: Comparison of profiles of reconstructed Shepp Logan image using different methods with 128 projections and a high SNR=40dB, after 50 iterations. The top image is the profile of original and reconstructed Shepp Logan image, the bottom image is a zone of the top image.

Figure 14 :

 14 Figure 14: Comparison of profiles of reconstructed Shepp Logan image using different methods with 128 projections and a low SNR=20dB, after 50 iterations. The top image is the profile of original and reconstructed Shepp Logan image, the bottom image is a zone of the top image.

Figure 15 :

 15 Figure 15: Comparison of the relative error of the reconstruction for the 64 × 64 Shepp Logan image (a) by using different methods with different projection numbers and a high SNR=40dB, and (b) by using different methods with different SNR and 64 projections. The relative error is measured after 50 iterations.

Figure 16 :

 16 Figure 16: Comparison of the relative error of the reconstruction for the 256 × 256 Shepp Logan image (a) by using different methods with different projection numbers and a high SNR=40dB, and (b) by using different methods with different SNR and 128 projections. The relative error is measured after 50 iterations.

Figure 17 :

 17 Figure 17: The final relative error by setting different value for λ. The top two figures show the QR method while the bottom two figures show the TV method. The left two figures present a high SNR case and the right two figures present a low SNR case.

Table 1 :

 1 1) and Table[START_REF] Arridge | Optical tomography in medical imaging[END_REF], different evaluations of the results are given, corresponding to the reconstruction results of Shepp Logan phantom of size 64 × 64 and 256 × 256 respectively. The evaluations are estimated with the reconstructed image within 50 iterations, using data with SN R = 20dB or SN R = 40dB. Different measurement of performance. The evaluations are estimated with 50 iterations for each method and each different case, for the Shepp Logan image of size 64 × 64. λ = 20 in QR and TV methods when there is a high SNR=40dB noise. λ = 300 in QR and λ = 40 in TV when there is a low SNR=20dB noise.

				64 projection					32 projection
		Data SNR=40dB	Data SNR=20dB	Data SNR=40dB	Data SNR=20dB
		QR	TV	HHBM	QR	TV	HHBM	QR	TV	HHBM	QR	TV	HHBM
	δ f	0.1164	0.1046	0.0255	0.1324	0.1146	0.1130	0.1812	1675	0.1114	0.1958	0.1761	0.2262
	ISNR	2.1171	2.5803	8.7153	1.7312	2.3562	2.4175	0.6862	1.0266	2.8000	0.6252	1.0843	-0.0015
	PSNR	21.4654	21.9286	28.0635	20.9068	21.5319	21.5931	19.5432	19.8836	21.6570	19.2071	19.6662	18.5804

Table 2 :

 2 Different measurement of performance. The evaluations are estimated with 50 iterations for each method and each different case, for the Shepp Logan image of size 256 × 256. λ = 20 in QR and TV methods when there is a high SNR=40dB noise. λ = 300 in QR and λ = 40 in TV when there is a low SNR=20dB noise.

				128 projection					64 projection
		Data SNR=40dB	Data SNR=20dB	Data SNR=40dB	Data SNR=20dB
		QR	TV	HHBM	QR	TV	HHBM	QR	TV	HHBM	QR	TV	HHBM
	δ f	0.0318	0.0147	0.0120	0.0839	0.0649	0.0715	0.0666	0.0443	0.0376	0.1323	0.0996	0.1110
	ISNR	3.0647	6.4358	7.3075	1.3952	2.5113	2.0895	1.9484	3.7192	4.4326	1.7058	2.9398	2.4657
	PSNR	27.1421	30.5133	31.3850	22.9369	24.0530	31.3850	23.9383	25.7091	26.4225	20.9580	22.1920	21.7179

Table 3 :

 3 Different measurement of performance. The evaluations are estimated until the convergence criterion is reached, for the Shepp Logan image of size 256 × 256. λ = 20 in QR and TV methods when there is a high SNR=40dB noise. λ = 300 in QR and λ = 40 in TV when there is a low SNR=20dB noise.

				128 projections
		Data SNR=40dB	Data SNR=20dB
		QR	TV	HHBM	QR	TV	HHBM
	δ f	0.0386	0.0208	0.0174	0.0932	0.0889	0.0766
	ISNR	2.2345	4.9143	5.6945	0.9380	1.1393	1.7866
	PSNR	26.3199	28.9918	29.7719	22.4797	22.6811	23.3283
	SSIM	0.9996	0.9998	0.9998	0.9991	0.9992	0.9993
				64 projections
		Data SNR=40dB	Data SNR=20dB
		QR	TV	HHBM	QR	TV	HHBM
	δ f	0.0729	0.0472	0.0387	0.1281	0.1057	0.0946
	ISNR	0.5528	3.4429	4.3076	1.8463	2.6800	3.1635
	PSNR	23.5427	25.4328	26.2975	21.0985	21.9321	22.4157
	SSIM	0.9993	0.9996	0.9996	0.9988	0.9990	0.9991
				32 projections
		Data SNR=40dB	Data SNR=20dB
		QR	TV	HHBM	QR	TV	HHBM
	δ f	0.1521	0.1290	0.0784	0.1925	0.1660	0.1254
	ISNR	2.0943	2.8105	4.9711	2.5575	3.1988	4.4188
	PSNR	20.3511	21.0673	23.2280	19.3292	19.9706	21.1905
	SSIM	0.9986	0.9988	0.9993	0.9982	0.9985	0.9989

L. Wang, A. Mohammad-Djafari, N. Gac / X-ray CT image reconstruction using Haar transform-based Bayesian model

L. Wang, A. Mohammad-Djafari, N. Gac / X-ray CT image reconstruction using Haar transform-based Bayesian model

L. Wang, A. Mohammad-Djafari, N. Gac / X-ray CT image reconstruction using Haar transform-based Bayesian model

* Thanks to China Scholarship Council for funding

A.2. Calculation of parameter γ (k) z

The gradient descent step length for updating z is calculated by minimizing the criterion J z (z -γ z dJ z ):

where J z (z -γ z dJ z ) is simplified as:

In Eq.( 44) we have:

So that:

The gradient of the criterion is:

By setting ∂J z (z -γ z dJ z ) ∂γ z | γz = 0, we get: