Paolo Irrera

Isabelle Bloch

Maurice Delplanque
email: -mdelplanque@eos-imaging.com

A DENOISING METHOD FOR WHOLE-BODY LOW-DOSE X-RAY IMAGES WITH ADAPTABLE PARAMETER CONTROL

Keywords: Low dose X-Ray images, Poisson noise, Non Local Means filter, Adaptive parameters

A denoising method is proposed for full body X-ray images, acquired under low dose conditions. The suggested algorithm is based on a non local means filter adapted to the statistics of Poisson noise. A new feature of the method is to locally set the filtering parameters in order to denoise while preserving details in low absorption regions. Thus, we propose to adapt the denoising parameters for each pixel by exploiting a global noise level measure and the standard deviation image of the gradient magnitude. Quantitative and visual results on phantom and real images show the interest of the method, achieving the objectives.

INTRODUCTION

In this paper we propose a method to reduce the amount of noise present in full body X-ray medical images, acquired under low dose conditions. The statistical model of the noise is a mixture of Gaussian and Poisson distributions. However, when the number of photons collected by the sensor is low, the quantum component is predominant [START_REF] Luisier | Image Denoising in Mixed Poisson-Gaussian Noise[END_REF]. This occurs typically for acquisitions performed on obese patients or, more generally, for low-dose image acquisitions.

The method is applied to EOS images. EOS imaging is a French company based in Paris that uses Charpak's high energy particle detector to enable low dose X-ray acquisitions [START_REF] Dubousset | Skeletal and spinal imaging with EOS system[END_REF]. A full body image with no vertical distortion can be obtained through linear vertical scan of the detector in just one acquisition and without the need of stitching algorithms. The main application field of interest is orthopedics. Benchmark studies have shown that EOS system reduces the entrance dose for the patient compared to traditional digital radiography systems (CR, DR) with significant better image quality [START_REF] Deschênes | Diagnostic imaging of spinal deformities: reducing patients radiation dose with a new slot-scanning X-ray imager[END_REF]. Nevertheless, increasing patient thickness or decreasing dose increases the level of noise in the images. This is why we propose a method to denoise the images before applying other post-processing algorithms, such as enhancement in low contrast regions. Obviously the denoising process should avoid losing useful information, and a robust denoising technique has then to be chosen, with parameters adapted to our context.

Several Poisson denoising algorithms are based on the Anscombe variance-stabilizing transform (VST). This operation allows associating the original data with a Gaussian statistical model. It has been formulated by Anscombe in [START_REF] Anscombe | The transformation of Poisson binomial and negative-binomial data[END_REF] and later generalized by Murtagh et al. in [START_REF] Murtagh | Image restoration with noise suppression using a multiresolution support[END_REF] to adapt it to contexts where Gaussian and Poisson noises are mixed up. The final goal of this kind of method is to be able to exploit

The PhD work of P. Irrera is supported by a grant from ANRT.

well known Gaussian denoising methods. An example is given in [START_REF] Zhang | Wavelets, ridgelets, and curvelets for Poisson noise removal[END_REF] where multi-scale image representation is combined with the VST. The main drawback lies in the non-linear distortion of noise-free elements introduced by the transformation. Besides, a more fundamental limit is the non-existence of a perfect variance stabilizing transform for some distributions, like for example the Poisson one [START_REF] Deledalle | How to compare noisy patches? Patch similarity beyond Gaussian noise[END_REF]. For this reason the VST step should be avoided in our context.

Luiser et al. proposed in [START_REF] Luisier | Fast Interscale Wavelet Denoising of Poisson-Corrupted Images[END_REF] a wavelets-based method to denoise elements affected by Poisson noise, based on a powerful data-driven technique that exploits an unbiased estimate of the mean square error, called Poisson Unbiased Risk Estimate (PURE). This work was then improved in [START_REF] Luisier | Image Denoising in Mixed Poisson-Gaussian Noise[END_REF] generalizing to the case of Poisson random vector corrupted by Gaussian noise. Another approach, directly processing the original data set elements, was proposed by Deledalle et al. in [START_REF] Deledalle | How to compare noisy patches? Patch similarity beyond Gaussian noise[END_REF] and [START_REF] Deledalle | Poisson NL Means: Unsupervised Non Local Means for Poisson Noise[END_REF]. The main difference is that this approach is based on Non-Local (NL) means, a nonlinear spatial filter that is gaining an increasing popularity in medical image processing for its performance [START_REF] Vega | Efficient and robust nonlocal means denoising of MR data based on salient feature matching[END_REF].

In this paper, we propose a denoising method based on Poisson NL means filter convenient to process EOS images and, theoretically, any full body X-Ray image. The objective is to denoise the image while preserving details in low absorption regions, i.e. slightly affected by noise. In particular, the problem of preserving the bone textures while reducing the global amount of noise is addressed in this paper. This goal is achieved by designing an adaptive local setting of parameters, as an original feature of the proposed approach. The method is part of a larger project which focuses on the design of specific algorithms for low dose X-Ray images. The dose is a central matter in nowadays society as shown by an important medical literature on exposure and related risks of cancer, see for example [START_REF] Ronckers | Multiple diagnostic X-rays for spine deformities and risk of breast cancer[END_REF] and [START_REF] Smith-Blindman | Use of Diagnostic Imaging Studies and Associated Radiation Exposure for Patients Enrolled in Large Integrated Health Care Systems, 1996-2010[END_REF]. EOS imaging answers the demand of low dose acquisition thanks to its unique system, but at the same time image processing needs to be exploited in order to provide the best image quality and enable specific applications. The way the method is adapted to our domain of interest is presented in Section 2. In Section 3 results on phantom and real images are provided, along with visual and quantitative evaluations. Conclusions and future developments are outlined in Section 4.

POISSON NL-MEANS WITH ADAPTIVE PARAMETER CONTROL

Poisson NL means

The key idea of NL means is that local is not always synonym of similar. In this type of filter, the pixel values are combined according to an affinity or similarity measure and not to a simple neighborhood constraint. Thus, in a window Ω centered at any pixel x i , weights w i, j are defined for each pixel x j in Ω. Then, the denoised value λi of pixel x i is estimated using the following weighted mean of the noisy pixel values k(x j ):

λi = ∑ |Ω| j=1 w i, j k(x j ) ∑ |Ω| j=1 w i, j (1)
For robustness reasons [START_REF] Deledalle | Poisson NL Means: Unsupervised Non Local Means for Poisson Noise[END_REF], the weight w i j is based on a pixel similarity assessed between patches Ψ around the pixels x i and x j , and should be defined according to a suitable statistical model. Given two observations of a Poisson process, k 1 and k 2 , two hypotheses are considered: either they share the same underlying noise-free value, i.e. λ 1 = λ = λ 2 , or they do not, i.e. λ 1 = λ 2 . Thus, in the case of a Poisson distribution, the similarity function evaluates the likelihood ratio [START_REF] Kay | Fundamentals of statistical signal processing[END_REF] corresponding to the two previous hypotheses [START_REF] Deledalle | Poisson NL Means: Unsupervised Non Local Means for Poisson Noise[END_REF]:

f (k 1 , k 2 ) = -log max λ p(k 1 |λ 1 = λ )p(k 2 |λ 2 = λ ) max λ p(k 1 |λ 1 = λ ) max λ p(k 2 |λ 2 = λ ) = k 1 log k 1 + k 2 log k 2 -(k 1 + k 2 ) log k 1 + k 2 2
(2)

Note that this formula replaces the Euclidean distance used for Gaussian distributions. Finally, the similarity function in Equation 2, is the input of an exponential kernel of bandwidth α in order to obtain the weights w i, j ∈ [0, 1]:

w i, j = exp - ∑ x ′ i ∈Ψ i ∑ x j ′ ∈Ψ j f (k(x i ′ ), k(x j ′ )) α (3)
where Ψ i (respectively Ψ j ) denotes the patch centered at x i (respectively x j ).

Alpha Map

It has been shown that in the case of images with a low signal to noise ratio, the quality of the denoising can be improved using a smoothed version of the input [START_REF] Kervrann | Optimal spatial adaptation for patch-based image denoising[END_REF]. Hence in [START_REF] Deledalle | Poisson NL Means: Unsupervised Non Local Means for Poisson Noise[END_REF] Equation 3 is redefined using a pre-estimate of the noise-free image and includes an additional similarity term computed between the smoothed image values. This requires an additional filter parameter β , analogous to α. The main aspect of the method proposed in [START_REF] Deledalle | Poisson NL Means: Unsupervised Non Local Means for Poisson Noise[END_REF] is to set jointly the parameters that control the shape of the exponential, i.e. α and β . In order to achieve this task the authors estimate the two parameters following a risk minimization principle. In particular, PURE [START_REF] Luisier | Fast Interscale Wavelet Denoising of Poisson-Corrupted Images[END_REF] is used, since the mean square error cannot be used when no prior knowledge is available.

The proposed method differs from Deledalle's one [START_REF] Deledalle | Poisson NL Means: Unsupervised Non Local Means for Poisson Noise[END_REF] in two main aspects: (i) the smoothed version of the original image is not taken into account, and (ii) we create a map of local alpha values rather than exploiting a global value. The first decision is related to the fact that, as we noticed after a few tests, with EOS images the advantage obtained in performance is not that significant to justify the increasing computation time. As for the second aspect, in medical images it is not trivial to define a global denoising parameter like a threshold or α in our case. Indeed, without prior models, extremely localized structures could fade away along with noise in the filtering processing. Alike, these structures could be saved but with a disappointing effect in other noisy regions. The interest of locally adapted processing has been also highlighted in different contexts. For instance, in [START_REF] Charnigo | A Semi-Local Paradigm for Wavelet Denoising[END_REF] Charnigo et al. proposed a semi-local paradigm for wavelet denoising of PET images exploiting a division of the image into suitable blocks that are, then, individually denoised.

Here we want to define a matrix of values of alpha, i.e. the alpha map. We propose to exploit measures computed from the standard deviation image of the amplitude of the gradient (from here on referred to as std image σ ). The alpha map is defined as follows:

α

(x i ) = M(ηζ (σ (x i ))) ( 4
)
where M is an average filter, η is the noise level that gives a hint on the amount of noise present in the image, and ζ (.) is a curve function of the std image value at pixel coordinates x i . The noise level η is computed as the maximum position of the histogram of the std image. This is, basically, the standard deviation computed on the background and provides an indication on the global amount of noise. Thus, for a given image, η is a constant value that modulates ζ (.). This curve is composed of two growing branches that are interpolated using Bsplines to guarantee C 1 -regularity. The first part is convex and grows slower than the second one which is concave (see Figure 1). The shape of ζ (.) is controlled by four parameters. Three parameters define the initial, final and switch amplitudes. The fourth parameter is a percentage of the global energy of the std image computed on its histogram. This allows defining the end of the first branch. The shape of the ζ curve allows controlling the strength of denoising on the image. The first branch points to regions that correspond to the background and low absorption parts of the body. In this way, the structured details are preserved thanks to a low value of α that practically limits the action of the NL-means filter. Then in the second branch, α grows up quickly to exploit the Poisson NL-means and, thus, to provide accurate denoising. This non-linear mapping of the denoising parameter allows solving the problem of poor performance of global methods on textured images as already highlighted by Luiser et al. in [START_REF] Luisier | Fast Interscale Wavelet Denoising of Poisson-Corrupted Images[END_REF]. The alpha map is then post-processed via an average filter in order to avoid strong transitions and visible changes of spatial resolution in the output image. An example of alpha map, computed on the image displayed in Figure 2(a), is presented in Figure 2(b). For instance, this case shows clearly how our non-linear mapping of α performs well on a clinical image: low values (cold colors) correspond mainly to background and bone textures, while the highest values (hot colors) are on the pelvis that is a not very much structured zone affected by a significant amount of noise.

The output of the denoising algorithm may still be poor in details, so we have to enhance the contrast. We use in this case a proprietary algorithm based on decomposition in Laplace multiscales. Although denoising and contrast enhancement are antagonist processing steps, what we propose here is useful and relevant. Indeed, the contrast enhancement algorithm will benefit of a proper representation of the input image.

EXPERIMENTAL RESULTS

The main goal of this section is to evaluate the interest and usefulness of a denoising step before the contrast enhancement. Thus, we will compare the results obtained in a reconstruction pipeline 1 without any denoising, and results obtained with the proposed technique. The contrast enhancement step is the same in both cases. Two types of images are considered for our tests: a phantom (standard PHD5000 fluoro phantom) used for image quality evaluations (Figure 3) and patient images (as the one in Figure 2(a)). Besides, for the diagnostic case analyzed in detail, we compared our method with the PURE-LET denoising [START_REF] Luisier | Fast Interscale Wavelet Denoising of Poisson-Corrupted Images[END_REF] and Poisson NL-means with global denoising parameters [START_REF] Deledalle | Poisson NL Means: Unsupervised Non Local Means for Poisson Noise[END_REF]. Note that for the tests here conducted with the global Poisson NL-means, the pre-estimate noise free image used in [START_REF] Deledalle | Poisson NL Means: Unsupervised Non Local Means for Poisson Noise[END_REF] is neglected and we set α = η. The standard PHD5000 phantom is covered by a layer of 20cm of PMMA blocks in order to simulate the thickness of an adult human body. The peak kilo voltage of the X-Ray generator is fixed at 63 kV. Two parameters can be tuned: the X-Ray tube current and the exposure time. For five different settings of these parameters, the corresponding dose values were computed: 31.68 µG, 63.48 µG, 97.98 1 Image processing chain that processes the data acquired by the sensor in order to provide a display image that can be clinically interpreted. µG, 136.34 µG, and 306.35 µG. As for the denoising parameters, the window size is fixed to 5 × 7 and the patch size to 3 × 3. The shape of the ζ (.) curve, i.e. the alpha map (see Equation 4) is defined according to the acquisition parameters: at lower dose amounts, the noise is more important and a stronger denoising is necessary. Three quantitative measures are considered: Signal to Noise Ratio (SNR), Contrast to Noise Ratio (CNR) and Dynamic (DYN). The SNR is computed as the ratio of the mean intensity value in a supposedly constant region on its standard deviation (Figure 3 -zone 1). The CNR is estimated in the fifth grey disk of the MediTest phantom (Figure 3 -zone 2). Finally, the DYN is computed as the percentage (on the 65536 possible gray levels) of the difference between the mean intensities of zones 3 and 4 (Figure 3). A comparison between the results obtained using the classical and the proposed methods is reported in Figure 4. Using the proposed approach, the SNR gains a factor two on a linear scale. The dynamics is reduced by 6% between the highest and the lowest dose values rather than 11% using the classical method. The most remarkable improvement concerns the CNR: the disk that is identifiable with the original approach only at 306.35 µG is now visible all along the considered range of doses since the CNR is larger than one.

Let us now illustrate the method on a lower limbs image. We defined two regions of interest (ROI): the head of the femur (Figure 5) and the knee (Figure 6). For instance, the first ROI is important to identify the anterior and posterior acetabulum2 lips. Instead, the second ROI was chosen in order to verify that the method preserves the bone trabecular structures. The quality of the outputs on the pelvis ROI are quantitatively compared in Table 1. In addition to the SNR and CNR, to establish whether the loss of resolution is significant or not, we relied on a no reference sharpness measure based on Local Phase Coherence (LPC) proposed by Hassen et al. in [16]. This measure is bounded in the interval [0, 1] and we refer to it as IQ lpc . Both the PURE-LET and the global NL-means denoising provide very similar results to those obtained using the classical pipeline, i.e. the global parameters do not allow a sufficiently strong denoising in this ROI. On the contrary, the measures in Table 1 and the picture in Figure 5(b) show a significant improvement of the image quality thanks to the proposed method. The sharpness of the result obtained with the proposed method is reduced but not significantly at all. As for the knee region, the results and the relative sharpness measures (the only measure that really matters in this case) are reported in Figure 6. The global Poisson NL-means, in this case, is not good at all because the bone textures fade completely away offering a not realistic outcome. PURE-LET performs better for this ROI, which makes us guess that its global automatic parametrization [START_REF] Luisier | Fast Interscale Wavelet Denoising of Poisson-Corrupted Images[END_REF] tends, for our images, to be quite conservative. Anyway, the proposed method is still the better among the studied denoising techniques thanks to its better ability to preserve fine details. Finally, we are able to show that our method is suitable to optimally adapt the strength of the denoising to different spots overcoming the drawbacks of other denoising methods. Indeed, the algorithm helps reducing the amount of noise in thick areas, e.g. the pelvis, while preserving the bone details in thinner zones, e.g. the knee. Similar results have been obtained on a series of 20 clinical images.

CONCLUSIONS

In this paper, we proposed a new denoising method, by adapting the NL means filter for Poisson noise model introduced in [START_REF] Deledalle | Poisson NL Means: Unsupervised Non Local Means for Poisson Noise[END_REF] to medical X-Ray images and in particular to those acquired with an EOS system. The main feature of the proposed approach is a local control of the parameters. This allowed us to overcome the drawbacks of other Poisson denoising techniques thanks to the capacity of reducing the amount of noise while preserving the bone textures in low X-rays absorption regions. Quantitative and visual results on both phantom and real images have shown the efficiency of the proposed method. As future work, other potentialities of the technique could be exploited. For example, a contrast enhancement algorithm could be designed in order to reach a better compatibility with the technique introduced in this paper. Anyway, this paper represents a first step towards our main goal: an optimized combination of hardware acquisition system and software processing to provide better and most useful image quality achievable for the lower dose delivered to the patient.

Fig. 1 :

 1 Fig. 1: ζ curve, as a function of the std image values. Initial value = 0.01; shift value = 0.1; final value = 10; shift position = 30%.

Fig. 2 :

 2 Fig. 2: (a) Input image. (b) Alpha Map: cold colors → low α values to preserve the details, hot colors → strong α values to reduce noise.

Fig. 3 :

 3 Fig. 3: Standard PHD5000 phantom and regions used for assessing the image quality.

Fig. 4 :

 4 Fig. 4: Image quality metrics on the phantom image: left SNR, center CNR, right DYN. Dotted lines: no specific processing for the noise; continuous lines: proposed NL-means filter with adaptive parameters.

Fig. 5 :

 5 Fig. 5: Pelvis ROI of a diagnostic image: (a) Pelvis without denoising, (b) Pelvis after denoising using the proposed method.

Fig. 6 :

 6 Fig. 6: Knee ROI of a diagnostic image: (a) Pelvis without denoising IQ lpc = 0.941, (b) PURE-LET denoising [8] IQ lpc = 0.921, (c) Global Poisson NL-Means [9] IQ lpc = 0.909, (d) Proposed method IQ lpc = 0.937.

Table 1 :

 1 Pelvis ROI. Quantitative comparison of classical EOS pipeline, PURE-LET[START_REF] Luisier | Fast Interscale Wavelet Denoising of Poisson-Corrupted Images[END_REF], global Poisson NL-means[START_REF] Deledalle | Poisson NL Means: Unsupervised Non Local Means for Poisson Noise[END_REF] and proposed method.

		Classical EOS PURE-LET [8] Global NL-means [9] Proposed
	SNR	18.55	19.87	20.52	21.19
	CNR	1.89	2.14	2.24	3.08
	IQ lpc	0.977	0.976	0.973	0.965
		(a)		(b)

Cup-shaped cavity of the pelvis into which the head of the femur fits


Jilles S Dibangoye

Olivier Buffet

Franc ¸ois Charpillet

Error-bounded Approximations for Infinite-Horizon Discounted Decentralized POMDPs

Keywords:

We address decentralized stochastic control problems represented as decentralized partially observable Markov decision processes (Dec-POMDPs). This formalism provides a general model for decision-making under uncertainty in cooperative, decentralized settings, but the worst-case complexity makes it difficult to solve optimally (NEXP-complete). Recent advances suggest recasting Dec-POMDPs into continuous-state and deterministic MDPs. In this form, however, states and actions are embedded into high-dimensional spaces, making accurate estimate of states and greedy selection of actions intractable for all but trivial-sized problems. The primary contribution of this paper is the first framework for error-monitoring during approximate estimation of states and selection of actions. Such a framework permits us to convert state-of-the-art exact methods into error-bounded algorithms, which results in a scalability increase as demonstrated by experiments over problems of unprecedented sizes.

Learning and planning algorithms for decentralized stochastic control problems are of importance in a number of practical domains such as network communications and control; rescue, surveillance and exploration tasks; multi-robotics; collaborative games [START_REF] Oliehoek | Heuristic search for identical payoff Bayesian games[END_REF]; to cite a few. Decentralized partially observable Markov decision processes (Dec-POMDPs) have emerged as a standard framework for modeling and solving such problems [START_REF] Bernstein | The complexity of decentralized control of Markov decision processes[END_REF]. This formalism involves a set of agents with different, but related, observations about the world, which cooperate to achieve a common long-term goal, but cannot explicitly communicate with one another. While many decentralized stochastic control problems can be formalized as Dec-POMDPs, only a few of them can be solved optimally due to their worst-case complexity: finite horizon problems are in NEXP, and infinite horizon problems are undecidable [START_REF] Bernstein | The complexity of decentralized control of Markov decision processes[END_REF]. This intractability is due to the doublyexponential growth in required computational resources, making it hard to find an optimal solution for all but the smallest instances [START_REF] Hansen | Dynamic programming for partially observable stochastic games[END_REF][START_REF] Bernstein | Policy iteration for decentralized control of Markov decision processes[END_REF].

A recent scalability increase builds upon two fundamental results [START_REF] Dibangoye | Optimally solving Dec-POMDPs as continuous-state MDPs[END_REF]. The first result establishes that Dec-POMDPs can be transformed with no loss of optimality into continuous-state and deterministic MDPs, called occupancy MDPs. In this form, the states -called occupancy states-are distributions over the states and actionobservation joint histories of the original Dec-POMDPs, and the actions -called decentralized decision rules-are mappings from joint histories to joint actions of the original Dec-POMDPs. Secondly, the optimal value function of a finite-horizon occupancy MDP is a piecewise linear and convex function of the occupancy state. These results allow to combine advances in continuous-state MDP and POMDP algorithms, which (among others) result in the feature-based heuristic search value iteration algorithm (FB-HSVI). This algorithm can produce optimal solutions for medium-sized problems and medium planning horizons, but quickly runs out of time and memory for larger-scale problems and planning horizons. Such limited scalability is mainly because states and actions of occupancy MDPs are embedded into high-dimensional spaces, making accurate estimate of states and greedy selection of actions intractable for all but trivial-sized problems.

A natural question to ask is whether approximate (error-bounded) solutions can be found efficiently for decentralized stochastic control problems. On the one hand, memory-bounded dynamic programming algorithms for solving infinite-horizon discounted Dec-POMDPs are often quite effective at finding good heuristic solutions, while requiring bounded computational resources [START_REF] Seuken | Formal models and algorithms for decentralized decision making under uncertainty[END_REF][START_REF] Dibangoye | Point-based incremental pruning heuristic for solving finite-horizon Dec-POMDPs[END_REF][START_REF] Kumar | Point-based backup for decentralized POMDPs: Complexity and new algorithms[END_REF][START_REF] Dibangoye | Toward error-bounded algorithms for infinite-horizon Dec-POMDPs[END_REF]. However, these methods do not come with rigorous guarantees concerning the quality of the final heuristic solution. On the other hand, error-bounded algorithms for solving infinite-horizon discounted Dec-POMDPs exist. Examples include error-bounded methods for discounted POMDPs that are (or can be) transferred back to discounted Dec-POMDPs: policy iteration (PI) [START_REF] Bernstein | Policy iteration for decentralized control of Markov decision processes[END_REF]; incremental policy generation (IPG) [START_REF] Amato | Incremental policy generation for finite-horizon Dec-POMDPs[END_REF]; point-based value iteration (PBVI) [START_REF] Macdermed | Point based value iteration with optimal belief compression for Dec-POMDPs[END_REF][START_REF] Pineau | Point-based value iteration: An anytime algorithm for POMDPs[END_REF]; and heuristic search value iteration (HSVI) [START_REF] Dibangoye | Optimally solving Dec-POMDPs as continuous-state MDPs[END_REF][START_REF] Smith | Point-based POMDP algorithms: Improved analysis and implementation[END_REF]. These algorithms rely either on ǫ-pruning methods 1 (PI and IPG) or/and on exploration strategies that focus on a small subset of the search space (PBVI and HSVI), but they all make use of greedy action-selection and accurate state-estimation operators, which quickly exhausts the available resources before convergence. Furthermore, theoretical analyses of point-based approaches (e.g., PBVI) demonstrate that resulting error bounds are loose and have only a theoretical significance [START_REF] Smith | Point-based POMDP algorithms: Improved analysis and implementation[END_REF][START_REF] Dibangoye | Toward error-bounded algorithms for infinite-horizon Dec-POMDPs[END_REF].

In this paper, we focus on characterizing efficient error-bounded solutions for infinitehorizon discounted decentralized stochastic control problems. The novel approach proceeds by converting infinite-horizon discounted Dec-POMDPs into finite-horizon discounted occupancy MDPs, thereby computing a non-stationary policy over a finite planning horizon. In such a setting, approximations are typically achieved by replacing greedy action-selection and accurate state-estimation operators by approximate counterparts. In addition, we preserve the ability to bound the error with respect to the optimal infinite-horizon value function. Our study differs from previous studies in that it directly bounds the regret, avoiding the max-norm machinery of previous analyses of value and policy iteration algorithms [START_REF] Puterman | Markov Decision Processes, Discrete Stochastic Dynamic Programming[END_REF][START_REF] Bagnell | Policy search by dynamic programming[END_REF][START_REF] Pineau | Point-based value iteration: An anytime algorithm for POMDPs[END_REF][START_REF] Smith | Point-based POMDP algorithms: Improved analysis and implementation[END_REF][START_REF] Dibangoye | Toward error-bounded algorithms for infinite-horizon Dec-POMDPs[END_REF][START_REF] Scherrer | Improved and generalized upper bounds on the complexity of policy iteration[END_REF], which may result in tighter bounds. We further extend the state-of-the-art feature-based heuristic search value iteration algorithm to incorporate the error we make in both greedy action-selection and accurate state-estimation. The result is an algorithm that can solve problems of unprecedented sizes from the literature while providing strong theoretical guarantees.

In the remainder of this paper, we will introduce in Section 1 the Dec-POMDP framework and the reformulation into occupancy MDPs. Section 2 extends, from finite-horizon settings to infinite-horizon ones, recent advances in optimally solving Dec-POMDPs as occupancy MDPs. Then, we describe the novel approximation framework, derive theoretical guarantees and algorithmic extensions in Section 3. Finally, we present in Section 4 experimental results demonstrating the scalability of the resulting error-bounded algorithm.

1 From Dec-POMDPs to OMDPs This section presents formalisms for the infinite-horizon discounted decentralized partially observable Markov decision process (Dec-POMDP) and its associated occupancy Markov decision process (OMDP).

Decentralized stochastic control problems as Dec-POMDPs

The Dec-POMDP framework formalizes a discrete stochastic system that evolves under the influence of N agents. A key assumption in this framework is that agents cannot directly observe the true state of the system. In fact, they have different but related observations about the state of the system and cannot explicitly communicate with one another. Nevertheless, they need to cooperate in order to achieve a common long-term goal, i.e., to select actions that maximize the collection of rewards in the long run.

Definition 1 (Dec-POMDP). An N-agent decentralized partially observable Markov decision process is given as a tuple M ≡ (S , {A i }, {Z i }, p, r, b 0 , γ), where: S is a finite set of hidden states; A i is a finite set of private actions of agent i ∈ {1, 2, . . . , N}; Z i is a finite set of private observations of agent i ∈ {1, 2, . . . , N}; p a,z (s, s ′ ) = Pr(s ′ , z|s, a) is a dynamics model of the team of agents as a whole; r(s, a) is a reward model of the team of agents as a whole; b 0 is an initial belief state; and γ ∈ (0, 1) is a discount factor.

The goal of solving M is to find an N-tuple of private policies π ≡ (π i ) i∈{1,2,...,N} that yields the highest discounted total reward starting at b 0 :

V π M,γ,0 (b 0 ) = E ∞ t=0 γ t • r(s t , a t ) | π, b 0 . (1)
Let decentralized policy π be a N-tuple of private policies (π i ) i∈{1,2,...,N} . Each private policy π i is a sequence of private decision rules (π i t ) t∈{0,1,...,∞} . The t-th private decision rule π i t : Θ i t → A i of agent i prescribes private actions based on the whole information available to the agent up to time step t, namely its complete history of past actions and observations

θ i t = (a i 0 , z i 1 , . . . , a i t-1 , z i t ) ∈ Θ i t .
We define Θ i t to be the set of all length-t private histories of actions and observations agent i may have experienced, Θ t ≡ × i∈{1,2,...,N} Θ i t the set of joint histories and Θ = ∪ t∈{0,1,...,T -1} Θ t . In addition, we define the t-th decentralized decision rule π t to be an N-tuple (π i t ) i∈{1,2,...,N} of private decision rules.

Since the history length grows as time goes on, for infinite horizon cases, this would require private decision rules to have infinite memory, which is not possible in practice. Therefore, we shall specify the nature of the decentralized policies we target in more detail. We first notice that the optimal value function over an infinite horizon can be arbitrarily accurately approximated by the optimal value function over a finite horizon.

To this end, we choose finite horizon T so that the regret of operating only over T = ⌈log γ ((1γ)ε/ r ∞ )⌉ steps instead of an infinite number of steps is upper-bounded by any arbitrarily small scalar ε > 0, where r ∞ = max{|r(s, a)| : ∀s ∈ S , ∀a ∈ A}. Indeed, the regret is upper-bounded by the cumulated sum of discounted losses from time step T onwards, so that: ∞ t=T γ t r ∞ ≤ ε. In the remainder of this paper, we restrict the search space to decentralized policies described over planning horizon T . Unlike infinite-horizon decentralized policies, finite-horizon decentralized policies require a finite memory. At the execution phase, agents follow actions their private policies prescribe up to time step T ; thereafter they behave randomly. Doing so, we are guaranteed to achieve performance with bounded error as discussed later below. Before proceeding any further, we next consider a reformulation of finite-horizon Dec-POMDPs into occupancy MDPs.

Occupancy Markov decision processes

The decentralized partially observable Markov decision process framework formalizes a decentralized stochastic control problem from a perspective oriented towards agents. In such a setting, agents are unaware of which actions the other agents take and which observations they receive; each agent behavior is based only upon its private histories. In this section, however, we formalize decentralized stochastic control problems from a perspective oriented towards centralized solution methods. In such a perspective, the system evolves under the control of agents based upon the total information about the state of the system the centralized solution method makes available to all agents prior to the execution phase, namely the information state.

The t-th information state ζ t ≡ (b 0 , π 0 , . . . , π t-1 ) is a sequence of decentralized decision rules starting at initial belief state b 0 . It satisfies the following recursion: ζ 0 ≡ (b 0 ) and ζ t ≡ (ζ t-1 , π t-1 ), for all t ∈ {1, . . . , T -1}. Next, it will prove useful to introduce the concept of occupancy states, as a means of maintaining a concise representation of the information state. A t-th occupancy state ξ t is a distribution P ζ t (s t , θ t ) over histories and hidden states of M conditional on an information state ζ t . For the sake of simplicity, we use notation Θ(ξ t ) to represent histories that are reachable in occupancy state ξ t . The occupancy state has many important properties. First, it is a sufficient statistic of the information state when estimating the (current and future) reward to be gained by executing a decentralized decision rule: R(ξ t , π t ) = s t θ t ξ t (s t , θ t ) • r(s t , π t (θ t )). In addition, it describes a deterministic and Markov decision process, where next occupancy state ξ t+1 ≡ P(ξ t , π t ) depends only upon the current occupancy state ξ t and the next decentralized decision rule π t :

ξ t+1 (s ′ , (θ t , a t , z t+1 )) = 1 {a t } (π t (θ t )) s∈S ξ t (s, θ t ) • p a t ,z t+1 (s, s ′ ), (2)
for s ′ ∈ S , a t ∈ A, z t+1 ∈ Z, θ t ∈ Θ and where 1 F is an indicator function. This process is known as the occupancy Markov decision process.

Definition 2 (OMDP). Let M ≡ (△, A, R, P, γ, ξ 0 , T ) be the T -steps OMDP with respect to Dec-POMDP M, where γ is a discount factor; ξ 0 corresponds to the initial belief in M; △ ≡ ∪ t∈{0,1,...,T } △ t is the set of occupancy states up to time T ; A ≡ ∪ t∈{0,1,...,T } A t is the finite set of decentralized decision rules; R(ξ t , π t ) is the reward model; P(ξ t , π t ) is the transition rule; and T is a planning horizon.

It is worth noticing that OMDP M can be seen as a generative model for occupancy states P(ξ t , π t ) and rewards R(ξ t , π t ), for all time step t ∈ {0, 1, . . . , T -1}. A recent result shows that an optimal solution for M, together with the correct estimation of the occupancy states, will give rise to the optimal solution of the original Dec-POMDP M over finite horizon T [START_REF] Dibangoye | Optimally solving Dec-POMDPs as continuous-state MDPs: Theory and algorithms[END_REF].

Optimally Solving Dec-POMDPs as OMDPs

This section reviews how to optimally solve Dec-POMDPs as OMDPs, a theory originally introduced under the total reward criterion [START_REF] Dibangoye | Optimally solving Dec-POMDPs as continuous-state MDPs[END_REF][START_REF] Dibangoye | Optimally solving Dec-POMDPs as continuous-state MDPs: Theory and algorithms[END_REF]. Here, we extend it to deal with the discounted total reward criterion.

Bellman's optimality equations

In this subsection, we extend dynamic programming properties, including Bellman's optimality equations, to OMDPs (respectively Dec-POMDPs). Before proceeding any further, we start with preliminary definitions.

The discounted total reward of a decentralized policy π ≡ (π t ) t∈{0,1,...,T -1} over T time steps and starting at occupancy state ξ t is

V π M,γ,t (ξ t ) = T -1 k=t γ k-t R(ξ k , π k ) | ξ k+1 = P(ξ k , π k ) , (3)
where the occupancy state sequence (ξ k ) k∈{t,t+1,...,T -1} is generated by the deterministic transition rule P under decentralized policy π: ξ k+1 = P(ξ k , π k ), ∀k ∈ {t, t + 1, . . . , T -1} and ∀t ∈ {0, 1, . . . , T -1}. Therefore, the optimal value function starting at occupancy state ξ 0 is V * M,γ,0 (ξ 0 ) = max π V π M,γ,0 (ξ 0 ). Hence, the optimal value function (V * M,γ,t ) t∈{0,1,...,T } is a solution of Bellman's optimality equation for M, given by:

V * M,γ,t (ξ t ) = max π t ∈A t R(ξ t , π t ) + γV * M,γ,t+1 (P(ξ t , π t )) , ∀ξ t ∈ △ (4)
and for t = T , we add a boundary condition V * M,γ,T (•) = 0. If it can be solved for (V * M,γ,t ) t∈{0,1,...,T } , an optimal decentralized policy π * ≡ (π * t ) t∈{0,1,...,T -1} may typically be obtained by maximization of the right-hand side for each ξ t , i.e.,

π * t ∈ arg max π t ∈A t R(ξ t , π t ) + γV * M,γ,t+1 (P(ξ t , π t )) , ∀ξ t ∈ △. (5)

Dynamic programming update operators

This subsection formally introduces the dynamic programming update operators involved in solving OMDPs, including: Bayesian state estimation; Bellman's evaluation and backup operators; and greedy action selection. To better understand this, let V be the set of real-valued functions f : △ t → R for all t ∈ {0, 1, . . . , T }.

Definition 3 (Bellman's evaluation operator). For each decentralized decision rule π t ∈ A t , let T π t : V → V be Bellman's evaluation operator, given by:

(T π t V M,γ,t+1 )(ξ t ) = R(ξ t , π t ) + γV M,γ,t+1 (P(ξ t , π t )), ∀ξ t ∈ △, π t ∈ A t . (6)
Bellman's evaluation operator transforms any arbitrary value function into a new value function based on a specified decentralized decision rule. It is worth noticing that Bellman's optimality equations (Equation 4) and greedy decision rule selections (Equation 5) can be stated in terms of the expression depending on occupancy state, decentralized decision rule and Bellman's evaluation operator. In the following, we formally define greedy selection and Bellman's update operators.

Definition 4 (Greedy action-selection operator). For each decentralized decision rule π t ∈ A t , let G : V → (△ → A) be the greedy operator, given by:

(GV M,γ,t+1 )(ξ t ) = arg max π t ∈A t (T π t V M,γ,t+1 )(ξ t ), ∀ξ t ∈ △, V M,γ,t+1 ∈ V. (7)
Together the greedy action-selection and Bellman's evaluation operators permit us to define Bellman's update operator as follows.

Definition 5 (Bellman's update operator). Let T : V → V be Bellman's update operator, given by:

(TV M,γ,t+1 )(ξ t ) = (T (GV M,γ,t+1 )(ξ t ) V M,γ,t+1 )(ξ t ), ∀ξ t ∈ △, V M,γ,t+1 ∈ V. (8)
Bellman's update operator maintains the value of a given occupancy state based on the greedy decentralized decision rule for a specified value function. When optimized exactly, the value function, solution of Bellman's optimality equations (Equation 4), is a piecewise-linear and convex function of the occupancy states [START_REF] Dibangoye | Optimally solving Dec-POMDPs as continuous-state MDPs: Theory and algorithms[END_REF]. That is, there exist finite sets of hyperplanes (Λ t ) t∈{0,1,...,T -1} , such that:

V * M,γ,t (ξ t ) = max λ t ∈Λ t s t ,θ t λ t (s t , θ t )• ξ t (s t , θ t ), where λ t ∈ R |S ||Θ t | for all t ∈ {0, 1, . . . , T -1}.
Mappings G and T serve to define a dynamic programming methodology for the solution of occupancy Markov decision process M. In particular, the piecewise-linearity and convexity property of the value function, together with mappings G and T, allow to combine advances in continuous-state MDP and POMDP algorithms, which have led to the development of a novel family of exact algorithms, including the feature-based heuristic search value iteration [START_REF] Dibangoye | Optimally solving Dec-POMDPs as continuous-state MDPs[END_REF][START_REF] Dibangoye | Optimally solving Dec-POMDPs as continuous-state MDPs: Theory and algorithms[END_REF].

The feature-based heuristic search value iteration

This subsection provides a succinct description of feature-based heuristic search value iteration (FB-HSVI) (Algorithm 1), which was originally introduced under the total reward criterion. Here, we extend it to address the discounted total reward criterion and discuss complexity issues.

The FB-HSVI algorithm's description. FB-HSVI extends to decentralized stochastic control problems the heuristic search value iteration (HSVI) algorithm, which was originally developed for partially observable Markov decision processes [START_REF] Smith | Point-based POMDP algorithms: Improved analysis and implementation[END_REF]. Similarly to HSVI, it corresponds to a family of trial-based algorithms that searches an optimal solution of an occupancy Markov decision process. FB-HSVI proceeds by generating trajectories of occupancy states, starting at the initial occupancy state. It maintains both upper and lower bounds over the optimal value function. It guides exploration towards occupancy states that are more relevant to the upper bound by greedily selecting decentralized decision rules with respect to the upper bound, and reducing the gap between bounds at visited occupancy states. If the gap between upper and lower bounds at the initial occupancy state is ε, then it terminates. In such a case, we are guaranteed FB-HSVI has converged to an ε-optimal solution, as initially targeted. Though FB-HSVI is already equipped with a mechanism for finding ε-optimal solutions -since it uses greedy action-selection and accurate state-estimation operators-in practice it quickly exhausts the available resources before convergence. To better understand this, we provide a complexity analysis of each operator involved in FB-HSVI.

Algorithm 1: The feature-based heuristic search value iteration for M (resp. M)

1 function FB-HSVI( M, ε, (V ¯M,γ,t ) t∈{0,1,••• ,T } , ( V M,γ,t ) t∈{0,1,••• ,T } ) 2 while Gap(ξ 0 ) > ε do Explore (ξ 0 ) 3 function Gap(ξ t ) 4 return V M,γ,t (ξ t )-V ¯M,γ,t (ξ t ) 5 function Explore (ξ t ) 6 if Gap(ξ t ) > ε/γ t then 7 π * t ← (G V M,γ,t+1 )(ξ t ) 8
Explore(P(ξ t , π * t ))

9

(T V M,γ,t+1 )(ξ t ) and (TV ¯M,γ,t+1 )(ξ t )

Complexity of dynamic programming operators. As FB-HSVI proceeds, there are three operations that can significantly affect the overall performance: the greedy actionselection operator G; the accurate state-estimation operator P; and finally, Bellman's update operator T. Given that the complexity of operators G, P and T are prohibitive for a number of realistic decentralized stochastic control problems, the importance of approximate variants is clear.

An error-bounded heuristic search framework

The primary contribution of this section is an error-bounded heuristic search framework which builds upon approximate variants of greedy action-selection and accurate state-estimation operators. We also provide a provable bound on the error FB-HSVI algorithms would make by using these approximate operators instead of their exact counterparts. The result is a general algorithmic framework that allows for monitoring the divergence between the exact and approximate solutions of infinite-horizon and discounted decentralized stochastic control problems represented as Dec-POMDPs.

Error-bounded action-selection operators

This subsection characterizes error-bounded action-selection operators that select decentralized decision rules within α of maximizing the value. Definition 6. Let α ∈ [0, ∞) T be a real vector. An α-approximate action-selection operator G : V → (△ → A) is such that, at each time step t ∈ {0, 1, . . . , T -1}, the decentralized decision rule found comes within α(t) of maximizing the value:

(T (GV M,γ,t+1 )(ξ t ) V M,γ,t+1 )(ξ t ) -(T ( GV M,γ,t+1 )(ξ t ) V M,γ,t+1 )(ξ t ) ≤ α(t), ∀ξ t ∈ △, V M,γ,t+1 ∈ V.
For any positive T -dimensional vector α, a feature-based heuristic search value iteration, together with an α-approximate action-selection operator, terminates with a final estimate V α M,γ,0 (ξ 0 ). The error between this approximate value and the optimal value is bounded and the bound depends only upon parameter α and γ.

Theorem 1. The error introduced in FB-HSVI by using G instead of G is bounded by

T -1 t=0 γ t α(t), assuming accurate estimation of the occupancy states during the planning phase. In particular, if α(t) = α(t + 1) = . . . = α for all t ∈ {0, 1, . . . , T -1}, then the error is bounded by 1-γ T 1-γ α.

Proof. Let π * and π α be decentralized policies that are optimal given that we use (P, G) and (P, G), respectively. Vectors ξ 1 , . . . , ξ T -1 being the occupancy states generated from ξ 0 when applying π * , it follows that:

V π * M,γ,0 (ξ 0 ) -V π α M,γ,0 (ξ 0 ) = T -1 t=0 γ t R(ξ t , π * t ) -V π α M,γ,0 (ξ 0 ) (definition of V π * M,γ,0 (ξ 0 )), = T -1 t=0 γ t R(ξ t , π * t ) + T -1 t=1 γ t V π α t:T -1 M,γ,t (ξ t ) -γ t V π α t:T -1 M,γ,t (ξ t ) -V π α M,γ,0 (ξ 0 ) (adding zero).
Next, we use the fact that V π α M,γ,T (ξ T ) = 0 to re-arrange terms:

= T -1 t=0 γ t R(ξ t , π * t ) + γ T V π α M,γ,T (ξ T ) + T -1 t=1 γ t V π α t:T -1 M,γ,t (ξ t ) -γ 0 V π α M,γ,0 (ξ 0 ) + T -1 t=1 γ t V π α t:T -1 M,γ,t (ξ t ) , = T -1 t=0 γ t R(ξ t , π * t ) + T -1 t=0 γ t+1 V π α t+1:T -1 M,γ,t+1 (P(ξ t , π * t )) -T -1 t=0 γ t V π α t:T -1 M,γ,t (ξ t ) , = T -1 t=0 γ t R(ξ t , π * t ) + γV π α t+1:T -1 M,γ,t+1 (P(ξ t , π * t )) -V π α t:T -1 M,γ,t (ξ t ) , = T -1 t=0 γ t V π * t ,π α t+1:T -1 M,γ,t (ξ t ) -V π α t:T -1 M,γ,t (ξ t ) , = T -1 t=0 γ t (T π * t V π α t+1:T -1 M,γ,t+1 )(ξ t ) -(T π α t V π α t+1:T -1 M,γ,t+1 )(ξ t ) , ≤ T -1 t=0 γ t (T (GV π α t+1:T -1 M,γ,t+1 )(ξ t ) V π α t+1:T -1 M,γ,t+1 )(ξ t ) -(T ( GV π α t+1:T -1 M,γ,t+1 )(ξ t ) V π α t+1:T -1 M,γ,t+1 )(ξ t ) , = T -1 t=0 γ t α(t).
Thus, the error between V π * M,γ,0 (ξ 0 ) and

V π α M,γ,0 (ξ 0 ) is bounded by T -1 t=0 γ t α(t). ⊓ ⊔
To the best of our knowledge, in decentralized stochastic control theory, this is the first attempt to monitor and bound the error made by using approximate action-selection instead of greedy action-selection. This bound comes with a natural interpretation: all time steps are not equally relevant to the final error. Indeed, due to discounted errors, approximate action-selection operators give more credit to errors they make at the earlier stages of the process. In other words, one can tolerate more approximation error at occupancy states that appear later in the process.

The problem of assigning errors to time steps goes beyond the scope of this paper, and will be addressed in the future. However, given the error vector α, another problem consists in finding a practical algorithm for selecting error-bounded actions over time steps. To do so, one can make use of the same branch-and-bound algorithms used for selecting greedy actions [START_REF] Dibangoye | Optimally solving Dec-POMDPs as continuous-state MDPs[END_REF][START_REF] Dibangoye | Optimally solving Dec-POMDPs as continuous-state MDPs: Theory and algorithms[END_REF]. Except that, now, these algorithms need to be interrupted whenever the gap between lower and upper bounds is α. In that case, we are guaranteed the returned action has value within α of the optimal value, as targeted.

Error-bounded state-estimation operators: definition and example

This subsection discusses the long term behavior of successive applications of an approximate state-estimation operator. Next, we formally define the family of approximate state-estimation operators we target. Then, we exhibit one such operator. And finally, we derive theoretical guarantees.

Since we are interested in quantifying the error between occupancy states, we choose the total variational distance as a metric for measuring their distance. The total variational distance between two probability distributions ξ and ξ

′ on [0, 1] |S ||Θ| is defined by ξ -ξ ′ TV = 1 2 s∈S ,θ∈Θ |ξ(s, θ) -ξ ′ (s, θ)|, ∀ξ, ξ ′ ∈ △.
Informally, the total variational distance ξξ ′ TV defines the minimal probability mass that would have to be re-assigned in order to transform occupancy state ξ into occupancy state ξ ′ . The following definition of approximate state-estimation operator P π t guarantees that, for any occupancy state ξ t ∈ △ t , we have ξ t P π tξ t P π t TV ≤ δ.

Definition 7. Let δ ∈ [0, 1] be a small scalar. Then, for each decentralized decision rule π t ∈ A t , transition matrix P π t is a δ-approximation of P π t if, for any occupancy state ξ t ∈ △ t , there exists δ ′ ∈ [0, δ] and (ξ ′ t+1 , ξ ′′ t+1 , ξ ′′ t+1 ) ∈ △ 3 t+1 such that

ξ t P π t = (1 -δ ′ )ξ ′ t+1 + δ ′ ξ ′′ t+1 and ξ t P π t = (1 -δ ′ )ξ ′ t+1 + δ ′ ξ ′′ t+1 .
Now, we introduce and describe Algorithm 2 for constructing an artificial occupancy state that is within δ (in terms of variational distance) from the original occupancy state. To ensure the total variational distance between artificial and original occupancy states is upper bounded by δ, the algorithm clusters together private histories of the original occupancy state that are close enough (see Definition 8). Then, it replaces each such cluster with a unique private history in that cluster. Finally, this private history represents the cluster in the artificial occupancy state.

Algorithm 2: The occupancy state approximation algorithm (OSA)

1 function OSA(ξ t , π t , δ) 2 ξ t+1 ← 0 and C ← Labels(ξ t , π t , δ)
C i ← ∅ and Θ i ← Θ i (ξ t P πt ) 8 while Θ i ∅ do 9 c i ← arg max θ i ∈Θ i |[θ i ] (ξt Pπ t ,δ) | 10 C i ← C i ∪ {c i } and Θ i ← Θ i \[c i ] (ξt Pπ t ,δ) 11 return ⊗ i∈{1,2,...,N} C i
Before proceeding any further, we introduce the criterion we use, namely the approximate probabilistic measure. Definition 8. Let ξ t be an occupancy state, and θ i and θi be two private histories in set Θ i (ξ t ). We say that θ i and θi are δ-probabilistically close if and only if:

Pr(X t , Y t |ξ t , θ i ) -Pr(X t , Y t |ξ t , θi ) TV ≤ δ, ( 9
)
where X t and Y t denote random variables associated with states and other agent histories, respectively. We also denote [θ i ] (ξ t ,δ) the entire set of private histories θi ∈ Θ(ξ t ) that are δ-probabilistically close to θ i and with respect to ξ t .

By clustering private histories that are δ-probabilistically close with a single private history in that cluster, we produce (from Definition 8) an approximate occupancy state ξt with respect to the original occupancy state ξ t such that: ξ t -ξt TV ≤ δ. Notice that Algorithm 2 is not guaranteed to produce an occupancy state with the minimum number of private histories. A more promising goal, which we do not address here, would be to find a clustering method that can identify the minimum number of clusters of private histories so that the total variational distance between original and artificial occupancy states is upper-bounded by δ.

Error-bounded state-estimation operators: theoretical analysis

We are now ready to bound the regret of using an approximate occupancy state instead of the accurate occupancy state. To do so, let P π 0:t-1 = P π 0 P π 1 • • • P π t-1 for all time steps t ∈ {1, 2, . . . , T }. Our analysis monitors the error we make step by step using approximate occupancy states.

Lemma 1. The total variational distance between ξ 0 P π 0:t-1 and ξ 0 P π 0:t-1 is bounded:

ξ 0 P π 0:t-1 -ξ 0 P π 0:t-1 TV ≤ 1 -(1 -δ) t , ∀t ∈ {1, 2, . . . , T }. ( 10
)
Proof. The proof holds directly by expanding ξ 0 P π 0:t and ξ 0 P π 0:t using Definition 7.

ξ 0 P π 0:t -ξ 0 P π 0:t TV , ≤ (1 -δ)ξ ′ 1 P π 1:t + δ ξ′′ 1 P π 1:t -(1 -δ)ξ ′ 1 P π 1:t -δξ ′′ 1 P π 1:t TV , ≤ (1 -δ) 2 ξ ′ 2 P π 2:t + δ(1 -δ) ξ′′ 2 P π 2:t + δ ξ′′ 1 P π 1:t -(1 -δ) 2 ξ ′ 2 P π 2:t -δ(1 -δ)ξ ′′ 2 P π 2:t -δξ ′′ 1 P π 1:t TV , ≤ (1 -δ) t ξ ′ t P π t:t + t k=1 δ(1 -δ) k-1 ξ′′ k P π k:t -(1 -δ) t ξ ′ t P π t:t -t k=1 δ(1 -δ) k-1 ξ ′′ k P π k:t TV , = t k=1 δ(1 -δ) k-1 ξ′′ k P π k:t -t k=1 δ(1 -δ) k-1 ξ ′′ k P π k:t TV , ≤ t k=1 δ(1 -δ) k-1 ξ′′ k P π k:t -ξ ′′ k P π k:t TV , ≤ t k=1 δ(1 -δ) k-1 , = 1 -(1 -δ) t . ⊓ ⊔
It is worth noticing that approximation errors tend to increase exponentially as time goes on. The following derives the regret induced by approximating state estimates. Theorem 2. Let δ ∈ [0, ∞) T be a scalar vector. The error introduced in FB-HSVI by using a δ-approximate state-estimation operator instead of the exact state-estimation operator is bounded by

2 r ∞ T -1 t=0 γ t [1 -t k=1 (1 -δ(k))],
assuming we use G for selecting decentralized decision rules. In particular, if δ(t) = δ for all time steps t ∈ {0, 1, . . . , T -1}, then the error is bounded by

2 1-γ T 1-γ -1-[γ(1-δ)] T 1-γ(1-δ) r ∞ .
Proof. Let π * and π be decentralized policies that are optimal given that we use accurate or approximate state-estimation operators, respectively. Define V π * M,γ,0 (ξ 0 ) as follows:

V π * M,γ,0 (ξ 0 ) = T -1 t=0 γ t R(ξ 0 P π * 0:t , π * t )
. Clearly, we have V π * M,γ,0 (ξ 0 ) ≤ V πM ,γ,0 (ξ 0 ) by definition of decentralized policy π. Using this property, we know that:

V π * M,γ,0 (ξ 0 ) -V πM ,γ,0 (ξ 0 ) ≤ V π * M,γ,0 (ξ 0 ) -V π * M,γ,0 (ξ 0 ), = T -1 t=0 γ t R(ξ 0 P π * 0:t-1 , π * t ) -R(ξ 0 P π * 0:t-1 , π * t ) .
Since the value function is piecewise-linear and convex, R π * t ≡ R(•, π * t ) is a linear function of occupancy states. Thus, if we let •, • be an inner product, then we have

V π * M,γ,0 (ξ 0 ) -V π * M,γ,0 (ξ 0 ) = T -1 t=0 γ t R π * t , ξ 0 P π * 0:t-1 -ξ 0 P π * 0:t-1 , (by linearity of R π * t ) ≤ T -1 t=0 γ t R π * t ∞ ξ 0 P π * 0:t-1 -ξ 0 P π * 0:t-1 1 , (Hölder's inequality) ≤ 2 T -1 t=0 γ t R π * t ∞ ξ 0 P π * 0:t-1 -ξ 0 P π * 0:t-1 TV , (where x 1 = 2 x TV ) ≤ 2 r ∞ T -1 t=0 γ t ξ 0 P π * 0:t-1 -ξ 0 P π * 0:t-1 TV , (where R π * t ∞ ≤ r ∞ ) ≤ 2 r ∞ T -1 t=0 γ t 1 -t k=1 (1 -δ(k)) .
This proves the result for any arbitrary δ ∈ [0, ∞) T . If we let δ(t) = δ for all time step t ∈ {0, 1, . . . , T -1}, then geometric series produce the following bound:

V π * M,γ,0 (ξ 0 ) -V πM ,γ,0 (ξ 0 ) ≤ 2 1-γ T 1-γ -1-[γ(1-δ)] T 1-γ(1-δ) r ∞ ,
which concludes the proof.

⊓ ⊔

Once again, in decentralized stochastic control settings, this is the first attempt to monitor and bound the error made by using approximate state-estimation operators. We note that, as time goes on, these operators become more tolerant to approximation errors. But there is no free lunch: approximation errors tend to increase as time goes on. This new bound provides a way to analyze this trade-off.

Convergence and error bounds

Given any arbitrary state-estimation and action-selection operators P and G, which come with provable guarantees, the feature-based heuristic search value iteration produces an estimate V M,γ,0 (ξ 0 ). The error between V M,γ,0 (ξ 0 ) and the true value function V * M,γ,0 (ξ 0 ) is bounded. The error depends on quantities ǫ, δ and α, each of which comes from a relaxation of the original problem. First, ǫ results from transforming an infinite horizon problem into a finite horizon one. Second, δ represents the vector of errors the state-estimation operator allows at each time step. Finally, α denotes the vector of errors the action-selection operator produces at each time step. Theorem 3. Let δ ∈ [0, ∞) T be the estimation operator parameter and α ∈ [0, ∞) T be the greedy operator parameter. The error of the feature-based heuristic search value iteration introduced by using P and G instead of P and G is bounded by:

2 r ∞ T -1 t=0 γ t        1 - t k=1 (1 -δ(k))        +         T -1 t=0 γ t α(t)         + ε, (11)
for any planning horizon T = ⌈log γ ((1γ)ε/ r ∞ )⌉.

Proof. Let π * , π α and π α,δ be decentralized policies that are optimal given that we use (P, G), (P, G) and ( P, G), respectively. Then,

V π * M,γ,0 (ξ 0 ) -V π α,δ M,γ,0 (ξ 0 ) = V π * M,γ,0 (ξ 0 ) -V π α M,γ,0 (ξ 0 ) + V π α M,γ,0 (ξ 0 ) -V π α,δ M,γ,0 (ξ 0 ) , ≤ 2 r ∞ T -1 t=0 γ t 1 -t k=1 (1 -δ(k)) + T -1 t=0 γ t α(t) . (12)
This bound together with the fact that we search only for T -step policies is sufficient to demonstrate that the result holds.

⊓ ⊔

This theorem provides the first result quantifying the influence of different approximate operators in the overall performance of an algorithm for solving Dec-POMDPs. To the best of our knowledge, no similar results exist in Dec-POMDPs.

Experiments

This section presents experiments on a selection of infinite-horizon γ-discounted Dec-POMDPs including small-sized benchmarks (broadcast channel, multi-agent tiger, recycling robots and meeting in a 3x3 grid) and large-sized benchmarks (box-pushing, mars rover and wireless). For each benchmark, we ran the error-bounded feature-based heuristic search value iteration (EB-FB-HSVI) algorithm using parameters ǫ (pruning criterion), α (action-selection tolerance), and δ (state-estimation tolerance). Notice that, over the selection of benchmarks, action-selection tolerance α has only minor influence on performance results, so we set α = 0 for many domains. We selected greedy actions using a constraint programming software, namely toulbar2 [START_REF] De Givry | Existential arc consistency: Getting closer to full arc consistency in weighted CSPs[END_REF]. EB-FB-HSVI ran on a Mac OSX machine with 2.4GHz Dual-Core Intel and 2GB of RAM available. 1. Results for infinite-horizon decentralized POMDPs with γ = 0.9, and by default we set ǫ = 0.001, α = 0 and δ = 0. Higher V ¯(ξ 0 ) is better. Results for Mealy NLP, EM, PeriEM, PI, MPBVI and IPG were likely computed on different platforms, an therefore time comparisons may be approximate at best. We compare EB-FB-HSVI for infinite-horizon Dec-POMDPs with state-of-the-art approximate and exact algorithms, including: optimal policy iteration (PI) [START_REF] Bernstein | Policy iteration for decentralized control of Markov decision processes[END_REF]; incremental policy iteration (IPG) [START_REF] Amato | Incremental policy generation for finite-horizon Dec-POMDPs[END_REF]; nonlinear programming (NLP and Mealy NLP) [START_REF] Amato | Optimizing fixed-size stochastic controllers for POMDPs and decentralized POMDPs[END_REF]; goal-directed algorithm [START_REF] Amato | Achieving goals in decentralized POMDPs[END_REF]; periodic expectation maximization algorithm (EM, Peri and PeriEM) [START_REF] Pajarinen | Periodic finite state controllers for efficient POMDP and DEC-POMDP planning[END_REF]; and modified point-based value iteration (MPBVI) [START_REF] Macdermed | Point based value iteration with optimal belief compression for Dec-POMDPs[END_REF]. Note that, while PI and IPG are optimal in theory, in practice they do not produce optimal solutions due to resources being exhausted before convergence. Table 1 reports performance results. For each domain and each algorithm, we report the lower-bound value function at the initial occupancy state V ¯(ξ 0 ), the computation time required to achieve that value, and the memory requirement |Λ|, which represents either the number of hyperplanes or the number of nodes in a policy graph.

Algorithm |Λ| Time V ¯(ξ 0 ) Broadcast (|S | = 4, |A i | = 2, |Z i | = 2) FB-HSVI 102
In all tested benchmarks, EB-FB-HSVI achieves values higher or equal to the highest values that have been recorded so far, while being multiple orders of magnitude faster than state-of-the-art algorithms over many domains. In particular, over smallsized problems, EB-FB-HSVI demonstrates the best trade-off between the quality of the solution and the computation time. In addition, it is the only algorithm to provide provable bounds on the resulting solutions. In the broadcast channel, for example, both EB-FB-HSVI and MPBVI provide the highest value known so far, but EB-FB-HSVI comes with two advantages over MPBVI. First, it guarantees that value 9.271 is within 0.001 of the optimum. Second, it computed this value four orders of magnitude faster than MPBVI. Over large-sized problems, EB-FB-HSVI terminated with the highest values over all benchmarks for parameters α = 0, δ = 0.01 and ǫ = 0.001. In the wireless problem, for example, the distance between the previous best value and EB-FB-HSVI's value is about 27, and EB-FB-HSVI was one order of magnitude faster than the previous best solver (MPBVI). 2. Theoretical guarantees of EB-FB-HSVI(δ = 0.01, ǫ = 0.001). We denote ε apriori the error computed a priori based on parameters δ and ǫ, and ε aposteriori the error computed a posteriori given approximation errors observed during the planning phase, both using Equation [START_REF] De Givry | Existential arc consistency: Getting closer to full arc consistency in weighted CSPs[END_REF]. Gap(ξ 0 ) = V ¯(ξ 0 ) -V(ξ 0 ), where V ¯(ξ 0 ) is provided by FB-HSVI(δ = 0.01, ǫ = 0.001) and V(ξ 0 ) results from EB-FB-HSVI(ǫ) for some ǫ.

We continue the study of the performance of EB-FB-HSVI with respect to tightness of error bounds. For each benchmark, we report in Table 2: a priori and a posteriori errors based on Equation [START_REF] De Givry | Existential arc consistency: Getting closer to full arc consistency in weighted CSPs[END_REF] for FB-HSVI(δ = 0.01, ǫ = 0.001); and gap Gap(ξ 0 ) =V ¯(ξ 0 ) -V(ξ 0 ) based on FB-HSVI(ǫ). Notice that a posteriori errors were computed based on approximation errors observed during the planning phase. Overall, a posteriori errors are tighter than a priori errors and closer to gaps. The tightness of a posteriori error is mainly because the observed approximation errors were significantly smaller than the targeted ones. In the tiger problem, for example, the a priori error is about 166.7 whereas the a posteriori error and the gap are close: 0.727 and 0.001, re-spectively. Surprisingly, in some domains such as mars rover and box-pushing, a posteriori errors are even smaller than gaps. This phenomenon occurs when EB-FB-HSVI(ǫ) exhausts the total available resources before convergence, i.e., the gap is larger than targeted error ǫ. The closeness between the gaps and the a posteriori errors demonstrate, at least over all tested domains, the tightness of our error bounds.

Discussion, conclusion and future work

This paper presented two relatively interdependent contributions towards error-bounded solutions for infinite-horizon discounted Dec-POMDPs. First, we introduce the first error-bounded algorithmic framework for monitoring and bounding the error we make by using approximate action-selection and state-estimation operators instead of their exact counterparts. Second, we extend the state-of-the-art algorithm for solving finitehorizon Dec-POMDPs, namely the feature-based heuristic search value iteration algorithm, to infinite-horizon discounted Dec-POMDPs. The major difference being that we can now use approximate operators instead of exact operators while still being able to provide theoretical guarantees on the quality of the resulting solution. Experimental results demonstrate that, when compared to state-of-the-art algorithms, the error-bounded feature-based heuristic search value iteration algorithm improves both values and computation times in many domains from the literature.

Though this paper provides the first attempts to monitor and bound the error made by using approximate operators in decentralized stochastic control, similar results exist in simpler settings. Such results can be traced back to max-norm-based analyses of value and policy iteration algorithms for γ-discounted MDPs [START_REF] Puterman | Markov Decision Processes, Discrete Stochastic Dynamic Programming[END_REF], which prove that for some error α at each iteration there exists a stationary policy within 2γ

(1-γ) 2 α of the optimum. This result led to the development of much research on convergence arguments for γ-discounted MDPs and extensions including partially observable cases [START_REF] Pineau | Point-based value iteration: An anytime algorithm for POMDPs[END_REF][START_REF] Smith | Point-based POMDP algorithms: Improved analysis and implementation[END_REF] and decentralized stochastic control settings [START_REF] Dibangoye | Toward error-bounded algorithms for infinite-horizon Dec-POMDPs[END_REF]. Closer to our performance guarantees, [START_REF] Scherrer | Improved and generalized upper bounds on the complexity of policy iteration[END_REF][START_REF] Bagnell | Policy search by dynamic programming[END_REF] developed variations of value and policy iteration algorithms for computing nonstationary policies in γ-discounted MDPs for which the performance bounds can be significantly improved by a factor of 1 1-γ . Hence, Theorem 1 can be viewed as an extension of [START_REF] Scherrer | Improved and generalized upper bounds on the complexity of policy iteration[END_REF] to decentralized stochastic control settings. However, Theorem 2 differs from previous performance bounds in many aspects. First, it is not derived from the max-norm analysis; instead we measure state-estimation errors we made steps by steps, which may result in tighter performance bounds. As a consequence, it does not fit within the standard scheme of performance bounds. Nonetheless, it allows us to accurately estimate errors made in practice on all tested benchmarks.

In the future, we plan to extend the feature-based heuristic search value iteration algorithm so as to learn how to dynamically assign approximation errors over time steps in order to minimize the total computation time while providing the targeted error bound. Another avenue we plan to follow, relies on how to automatically find the minimum number of clusters of private histories such that the artificial occupancy state based on clusters is within δ of the original occupancy state.

3 foreach 4 return ξ t+1 5 function

 345 s ∈ S and c ∈ C do ξ t+1 (s, c) ← θ∈[c] (ξt Pπ t ,δ) ξ t P πt (s, θ) Labels(ξ t , π t , δ) 6 foreach i ∈ {1, 2, . . . , N} do 7

  To better understand the complexity involved in these operations, let |V| be the size of value function V (respectively the upper-or lower-bound value functions). Let Θ i (ξ t ) be the set of private histories of agent i involved in occupancy stateξ t , |Θ * (ξ t )| = max i∈1,2,...,N |Θ i (ξ t )| and |A * | = max i∈1,2,...,N |A i |.Algorithm 1 (lines 7 and 9) performs a greedy action-selection operator G, which involves enumerating and evaluating exponentially many decentralized decision rules in the worst case, and requires time complexity O(|V| |Θ * (ξ t )| N|A * | ) that grows doubly exponentially with increasing number of private histories involved in the occupancy state ξ t . In practice, branch-and-bound methods explore only a small portion of this set, which saves considerable time[START_REF] Dibangoye | Optimally solving Dec-POMDPs as continuous-state MDPs[END_REF][START_REF] Dibangoye | Optimally solving Dec-POMDPs as continuous-state MDPs: Theory and algorithms[END_REF]. Then, Algorithm 1 (line 8) computes the next occupancy state given the current one and the next decentralized decision rule. Unlike the greedy action-selection operator, this state-estimation rule has complexity O(|S | 2 |Θ(ξ t )||Z|), that is polynomial in the number of joint histories involved in the occupancy states and the number of joint observations. However, in the worst case, the number of joint histories increases by a factor of |Z| as time goes on. This may limit ability to perform the greedy action-selection operator later on. Finally, Algorithm 1 (line 9) performs Bellman's update operator T to maintain both upper and lower bounds at a given occupancy state, namely point-based Bellman's update. Unlike the full Bellman's update operator, the point-based Bellman's update operator maintains the value function only at a single occupancy state at a time, which makes it significantly more tractable. Nonetheless, the complexity of this operation remains time demanding as it requires performing a greedy action selection.

  Recycling robots(|S | = 4, |A i | = 3, |Z i | = 2) Box-pushing (|S | = 100, |A i | = 4, |Z i | = 5)

		Algorithm	|Λ|	Time	V ¯(ξ 0 )
	19.8s 7.8s 36 < 18000s FB-HSVI(δ = 0.01) 435 MPBVI NLP 2 1s Dec-tiger (|S | = 2, |A i | = 3, |Z i | = 2) FB-HSVI(δ = 0.01) 52 6s 13.448 9.271 9.269 9.27 9.1 FB-HSVI 25 157.3s 13.448 MPBVI 231 < 18000s 13.448 Peri 10×30 220s 13.45 PeriEM 7×10 6540s 9.42 Goal-directed 11 75s 5.04 Mealy NLP 4 29s -1.49 EM 6 142s -16.3	FB-HSVI(δ = 0.01) FB-HSVI(α = 1, δ = 0.05) 288 331 FB-HSVI(ǫ = 30) 264 MPBVI 305 > 18000s 224.12 1715.1s 224.43 1405.7s 224.26 15.24s 199.42 Goal-directed 5 199s 149.85 Peri 15 × 30 5675s 148.65 Mealy NLP 4 774s 143.14 PeriEM 4 × 10 7164s 106.68 Mars rover (|S | = 256, |A i | = 6, |Z i | = 8) FB-HSVI(δ = 0.01) 136 74.31s 26.94 FB-HSVI(α = 0.2) 149 85.72s 26.92 FB-HSVI(ǫ = 1) 155 32.5s 26.77 Peri 10 × 30 6088s 24.13
	FB-HSVI FB-HSVI(δ = 0.01) 108 109 MPBVI 37 < 18000s 2.6s 0s Mealy NLP 1 0s Peri 6×30 77s PeriEM 6×10 272s EM 2 13s IPG 4759 5918s PI 15552 869s Meeting in a 3x3 grid (|S | = 81, |A i | = 5, |Z i | = 9) 31.929 31.928 31.929 31.928 31.84 31.80 31.50 28.10 27.20 FB-HSVI 108 67s 5.802 FB-HSVI(δ = 0.01) 88 45s 5.794 Peri 20×70 9714s 4.64	Goal-directed Mealy NLP PeriEM EM Wireless (|S | = 64, |A i | = 2, |Z i | = 6) 6 956s 3 396s 3 × 10 7132s 3 5096s FB-HSVI(δ = 0.01) 897 6309s -144.24 21.48 19.67 18.13 17.75 FB-HSVI(α = 0.1) 408 6740s -140.37 FB-HSVI(ǫ = 20) 866 6084s -176.59 MPBVI 374 > 18000s -167.10 EM 3 6886s -175.40 Peri 15 × 100 6492s -181.24 PeriEM 2 × 10 3557s -218.90 Mealy NLP 1 9s -294.50
	Table

An ǫ-pruning method circumvents regions of the search space that cannot significantly improve the current solution, and the resulting solution is guaranteed to be at ǫ of the optimum.


Paolo Irrera

Isabelle Bloch

Maurice Delplanque
email: -mdelplanque@eos-imaging.com

Contrast enhancement of Micro Dose X-ray images

Keywords: Micro Dose X-ray imaging, Multi-scale decomposition, Non-local means filter

teaching and research institutions in France or abroad, or from public or private research centers.

INTRODUCTION

Contrast enhancement is often applied on X-ray images because it improves the local contrast, thus easing the diagnosis. In radiology, while the large density variations occur mainly over low frequencies, the most relevant details belong to high frequencies [START_REF] Stahl | Noise-Resistant Weak-Structure Enhancement for Digital Radiography[END_REF]. The linear multi-scale (MS) methods, e.g. the Laplacian Pyramid [START_REF] Burt | The Laplacian Pyramid as a Compact Image Code[END_REF], allow separating an image in scales of spatial frequency. Thus, they can be used to enhance X-ray images [START_REF] Stahl | Noise-Resistant Weak-Structure Enhancement for Digital Radiography[END_REF], [START_REF] Dippel | Multiscale Contrast Enhancement for Radiographies: Laplacian Pyramid Versus Fast Wavelet Transform[END_REF]. However, they raise two issues. The first problem is the creation of halo artifacts due to smoothing of strong edges. In other image processing domains (e.g. digital photography), edge-preserving algorithms [START_REF] Paris | Local Laplacian Filters: Edge-aware Image Processing with a Laplacian Pyramid[END_REF] address this question. In radiology, these artifacts constitute an actual problem in presence of radio-opaque objects, e.g. prosthesis. Indeed, shadows may appear around the prosthesis making difficult for the radiologist to evaluate the health of the surrounding tissues. Excluding this case, the halo artifacts are very limited and can, thus, be tolerated [START_REF] Dippel | Multiscale Contrast Enhancement for Radiographies: Laplacian Pyramid Versus Fast Wavelet Transform[END_REF]. This is an assumption made in this paper. The second problem concerns the analysis of noisy images. The noise is represented at fine levels of the decomposition and, so, mixed up with details relevant for the diagnosis. Thus, the challenge consists in enhancing the fine details while avoiding boosting the noise. A common solution is the noise containment [START_REF] Stahl | Noise-Resistant Weak-Structure Enhancement for Digital Radiography[END_REF]: the noisy coefficients are not enhanced, but not attenuated either, which guarantees to preserve all the information.

In this project we process Micro Dose (MD) images acquired with the X-ray system EOS. This system is based on Charpak's high energy particle detector [START_REF] Dubousset | Skeletal and spinal imaging with EOS system[END_REF]. Frontal and lateral images of the full body of the patient can be acquired in upright position. EOS has been shown to provide spine radiographs at 6 to 9 times less dose than computer radiography [START_REF] Deschênes | Diagnostic imaging of spinal deformities: reducing patients radiation dose with a new slot-scanning X-ray imager[END_REF], which is why we call Low Dose (LD) the normally exposed EOS images. Micro Dose is a new EOS feature that allows a further dose reduction equal to 5 to 7 times w.r.t. LD acquisitions. The amount of dose delivered to the patient during a MD The PhD work of P. Irrera is supported by a grant from ANRT.

exam is comparable to 5 to 10 days of natural exposure. This feature is clearly of great interest given the increasing concerns on patients health due to over-exposed acquisitions [START_REF] Smith-Blindman | Use of Diagnostic Imaging Studies and Associated Radiation Exposure for Patients Enrolled in Large Integrated Health Care Systems, 1996-2010[END_REF]. These images are mainly useful for follow-up exams (e.g. measure of scoliosis angle trending). Nevertheless, the reduction of the flux of photons implies a rise of the amount of noise. So, since the noise propagates through coarser levels than the first one, the aforementioned noise containment technique is not optimal anymore. Then, a possible method consists in denoising the image before enhancing the local contrast [START_REF] Sakata | Noise Reduction and Contrast Enhancement for Small-Dose X-ray Images in Wavelet Domain[END_REF], [START_REF] Irrera | A Denoising Method for Whole-Body Low-Dose X-Ray Images with Adaptable Parameter Control[END_REF]. Nevertheless, this kind of approach could lead to an artificial look of structures at medium frequencies, e.g. the edges of the vertebrae. So, at present, MD images are processed in clinical routine with the same family of algorithms used on LD images.

Our contribution consists in optimizing the visualization of MD images in terms of: noise reduction, contrast enhancement and preservation of their natural aspect. We define a new MS decomposition that is not biased by the presence of noise that affects the input image. For this aim, we combine the decomposition of the original image and of its denoised version, including the denoising directly into the enhancement loop. The denoised image is obtained using a Non-Local means (NL-means) filter [START_REF] Buades | A Non-Local Algorithm for Image Denoising[END_REF], [START_REF] Deledalle | How to compare noisy patches? Patch similarity beyond Gaussian noise[END_REF] with an adaptable parameter control related to X-ray absorption that we described in [START_REF] Irrera | A Denoising Method for Whole-Body Low-Dose X-Ray Images with Adaptable Parameter Control[END_REF]. We want to stress that the proposed method is not just a MS denoising technique. Indeed, since we aim at obtaining a well contrasted output image, the reduction of the noise represents only half of the whole matter.

The method is explained in Section 2. In Section 3 we validate our algorithm by testing it on phantoms and MD images. Conclusions and future developments are outlined in Section 4.

CONTRAST ENHANCEMENT OF NOISY IMAGES

The algorithm is divided into three main steps (see Fig. 1) : the adaptive NL-means filter (block 1), the MS analysis (block 2) and the MS synthesis (block 3).

Adaptive NL-means denoising

The NL-means filter denoises an image through a representation in sub-regions called patches. The interest of working with patches lies in the robustness to noise and good representation of texture [START_REF] Deledalle | How to compare noisy patches? Patch similarity beyond Gaussian noise[END_REF]. Formally, given a noisy image I, the denoised pixel at the spatial location x i is estimated as a weighted sum:

Î (x i ) = ∑ |Ω| j=1 ς i, j I x j ∑ |Ω| j=1 ς i, j , ς i, j = exp - d P i , P j α ( 1
)
where Ω is the size of a window centered at pixel x i that determines which pixels have to be averaged. The weights ς i, j are defined in function of the similarity between the patches P i and P j centered respectively at the pixels x i and x j . The distance function d(.) is defined according to the type of noise that corrupts the image. It is generally the Euclidean distance, but other functions can be defined when the noise is not Gaussian [START_REF] Deledalle | How to compare noisy patches? Patch similarity beyond Gaussian noise[END_REF]. Finally, α is the bandwidth of the exponential kernel.

In a former work [START_REF] Irrera | A Denoising Method for Whole-Body Low-Dose X-Ray Images with Adaptable Parameter Control[END_REF], we pointed out that the choice of one value of α on a full body X-ray image is sub-optimal: there is a strong risk of loosing anatomical information in low absorption regions, e.g. bone textures. Thus, considering that the noise is stronger in low incident X-ray dose areas [START_REF] Stahl | Noise-Resistant Weak-Structure Enhancement for Digital Radiography[END_REF], we defined a map of different values empirically related to the thickness of the body of the patient. In this work we exploit the same idea. However the definition of the alpha map is simplified. First, the envelop of the body is segmented by thresholding the image in the linear domain 1 . Then, using the cumulative histogram of gray levels of the previously segmented image, we define a given number of regions of increasing thickness. Subsequently, we compute the local standard deviation image and label each region with the most frequent value in it, obtaining the noise level image. Finally, the alpha map is proportional to this noise level image.

Even if the noise model could be described as Poisson/Gaussian mixture, we employ an Euclidean distance to compare the patches. The signal dependent nature of the noise is taken into account through the alpha map that changes in function of the thickness. Thus, the estimation of the noise free image is practically unbiased.

Analysis of noise corrupted images

The analysis step is the core of the algorithm. The term analysis indicates the decomposition in scale levels S 0 I , S 1 I , . . . , S n I and detail levels

D 0 I , D 1 I , . . . , D n-1 I . The scale level S i I is the smoothed version of scale S i-1 I (with S 0 I = I), i.e. S i I = f i S i-1 I
where f i (.) is a Gaussian-like filter. We exploit an undecimated decomposition to avoid aliasing artifacts due to non-linear processing in the reconstruction. We use a simpler decomposition technique than the wavelets because this last one has been proven to be inappropriate for the contrast enhancement of medical images [START_REF] Dippel | Multiscale Contrast Enhancement for Radiographies: Laplacian Pyramid Versus Fast Wavelet Transform[END_REF]. The detail level is the output of a band-pass filter: D i-1 I = S i-1 I -S i I . We define the optimal decomposition using the levels obtained from the original image I and its denoised version Î. Precisely, the output detail level D i O is the weighted mean of the noisy detail level D i I and the denoised detail level D i Î . The decomposition of Î has the advantage of representing only details proper to the signal, avoiding the propagation of noise at levels i > 0. Nevertheless, the attenuation of coefficients can entail an excessive regularization of the outcome. However, thanks to the adaptive NL-means filter, we can preserve the information in low absorption regions and the edges are not smoothed. Then, we can detect these spots of the image and use the coefficients from the noisy decomposition to represent them. On the contrary, the flat regions are represented by the coefficients from the denoised image. This idea is formalized by the following weighting function:

w i = |∇(S I i )| |∇(S I i )| + |∇(S i Î )| β , β 1 (2)
where |∇(I)| is the gradient magnitude of an image I and β is a userdefined parameter. The value of w i is low when the gradient of the scale image obtained from the denoised image is higher than the gradient of the scale image obtained from the noisy image. This happens when an edge is preserved and the noise around it is reduced. In this case, we have detected an edge and, so, we give a higher weight to the coefficients from the noisy image. Conversely, high values of w i (.) imply a reduction of noise in smooth regions and the denoised coefficients are retained. The parameter β controls the power of the denoising: the higher it is, the lower are the weights. For example, assuming that the gradients are equivalent in a location x, if β = 1 then the coefficients have the same importance (w(x) = 0.5), instead if β = 2 the noisy coefficients are more relevant (w(x) = 0.25). Note that the values of w i (.) are lower and lower at coarser scales because the power of the noise gets less and less important. Finally, the output levels are obtained as

D i O = w i D i Î +(1-w i )D i I
, resulting in a denoised image but preserving the coefficients from strong edges.

Synthesis of noise corrupted images

The synthesis is the process that rebuilds the image from the detail levels and the low frequency (LF) residuals obtained during the analysis step. The enhanced image is obtained by summation:

I e = n-1 ∑ i=0 b i D i O + lut(S n Î ) (3)
The details D i O are boosted by the non linear function b i . In this paper, we use the boosting function defined by Stahl et al. [START_REF] Stahl | Noise-Resistant Weak-Structure Enhancement for Digital Radiography[END_REF]:

b i (d) = G i 1 - |d| d 0 p i + 1 if|d| d 0 1 if|d| > d 0 (4)
where G i is the maximum gain, p i controls how fast the amplification decreases toward 1 and d 0 is the amplitude value after which the coefficients are not enhanced. The value of d 0 is defined according to the binary levels used to encode the image. Instead, G i and p i are defined according to the scanned anatomical region and the application [START_REF] Stahl | Noise-Resistant Weak-Structure Enhancement for Digital Radiography[END_REF].

To treat the LF residual, we apply a sigmoid-like function (lut) that attenuates the low values of the LF residual and preserves the strong ones. We define γ ∈ [0, 1) as the strongest attenuation. The point where the function changes of convexity is fixed from the intersection of the cumulative histogram of the LF residual with the cumulative histogram of an ideal uniform distribution. On an X-ray image this allows maximizing the dynamic of the image.

VALIDATION OF THE METHOD

The proposed method was validated on phantoms and MD clinical images encoded on 16 bits. The images of the phantom portrait all the same object (standard PHD5000 fluoro phantom) but with different acquisition parameters that are: the peak kilo voltage output of the X-ray generator (kV), the X-ray tube current in mA (mA) and the speed factor of the C-scan (C∈ {1, 2, . . . , 8}, the higher it is the lower is the speed). The phantoms are used to quantify the improvement in terms of signal to noise ratio (SNR), contrast to noise ratio (CNR) and dynamic (DYN) that represents the global contrast. Our clinical data set includes 13 Micro Dose studies of spine acquired at the hospital Robert Debré, Paris 2 . The minimum and maximum values of dose (Air Kerma) are 7.37 µGy and 29.44 µGy and 17.43 µGy and 55.13 µGy for, respectively, frontal and lateral acquisitions. The parameters used in the algorithm are mostly defined in function of the content of the image. Indeed, the alpha map is proportional to the noise level map (see Sec 2.1); for all the results reported here the proportionality factor is equal to 4. The patches and the windows of the NL-means filter are respectively of size 7 × 7 and 21 × 21. In the analysis the only parameter to be defined is β (see Sec. 2.2). This parameter allows weighting the contribution of the denoising and it is the only one that needs to be changed as a function of dose. Thus, for our tests we used β = 1.5 and β = 3 respectively on MD and LD data. Finally, the synthesis parameters are practically defined on the basis of the type of study [START_REF] Stahl | Noise-Resistant Weak-Structure Enhancement for Digital Radiography[END_REF]. For all our tests, even those on phantoms, we opted for parameters that aim at enhancing the details of the spine. Thus, the high/medium frequency detail levels (i ∈ [START_REF] Stahl | Noise-Resistant Weak-Structure Enhancement for Digital Radiography[END_REF][START_REF] Dippel | Multiscale Contrast Enhancement for Radiographies: Laplacian Pyramid Versus Fast Wavelet Transform[END_REF]) are the most enhanced. Besides, we set γ = 0.5 to process the LF residual. The method was compared with two noise containment based algorithms that are usually used by EOS in clinical routine. One is employed on LD images and the other on MD images; they will be respectively named hereafter EOS-LD and EOS-MD. This last one is specifically designed to enhance the spine, which is why we chose this category of study.

Results on phantoms

The SNR is given by the ratio of the mean signal on the standard deviation in a constant region. The CNR is the ratio of the difference of signal between the first disk from the left in Fig. 2 and the background, on the standard deviation of this last one. Then, DYN is the difference of the average signal of the white and black tablets (see the upper zone of the image in Fig. 2). The results obtained with the proposed approach are better for all the metrics both in the LD and MD cases as shown, respectively, in Tab. [START_REF] Stahl | Noise-Resistant Weak-Structure Enhancement for Digital Radiography[END_REF] We would like to thank the Professor G. Sebag and the Doctor M. Alison from the hospital Robert Debré for providing the anonymized RAW images and agreeing to use them for our research purposes LD images, the SNR, the CNR and the DYN are improved on average by, respectively, ×1.80, ×3.72 and 8 points. Besides, the details are totally preserved as observed in Fig. 2, e.g. see the high contrast resolution grid and the writing. These results suggest that even in some cases with a sufficient amount of dose, the methods based on noise containment can fail in avoiding the noise boosting. Instead, the method that we propose proves to be more robust. As for the MD However, the scores of the DYN reveal a much better global contrast offered by our solution. Finally, the spatial resolution is not reduced in this case either.

Results on Micro Dose spine images

The validation of the method on clinical images was done first by a subjective comparison between the three enhancement algorithms.

In details, we study the lumbar region cropped from one of the MD images composing our clinical data set. The image enhanced with EOS-LD (Fig. 3(b)) has not an ideal global contrast. For example, the spine and the surrounding soft tissues have nearly the same gray level. Besides, the presence of noise is evident. This probably implies that the medium frequency details are not sufficiently enhanced and that the noise containment is not efficient. The solution offered by EOS-MD (Fig. 3(c)) is better, which justifies the current choice of this algorithm in the MD context. However, we can clearly see background noise in smooth regions. This is due to the inability of the MS decomposition to confine the noise at the finest level. Our algorithm (Fig. 3(d)) allows smoothing the most regular regions by reducing the visibility of the noise, while preserving the shapes of the vertebrae, i.e. only the relevant information of the image is enhanced. Finally, the global contrast is well balanced. We used Average Local Variance (ALV) metrics as objective assessment [START_REF] Chang | Image Contrast Enhancement Based on a Histogram Transformation of Local Standard Deviation[END_REF]. This approach consists in, first, classifying each pixel in three categories of regions, i.e. smooth, details and edges. Then, the average standard deviation is computed in each region obtaining, respectively, ALVS, ALVD and ALVE. Generally, the best contrast enhancement algorithm has low values of ALVS and ALVE, which quantify respectively the noise enhancement and the presence of halos, and high ALVD, which quantifies the details boosting. For every image of the clinical data set, we defined the three regions by a basic threshold classification of the detail levels D i Î and we computed the scores on the input image and on the outputs of the compared contrast enhancement algorithms. The averages of the metrics obtained on the whole data set are reported in Tab. 3. The algorithms EOS-LD and EOS-MD have similar scores: the significant enhancement of the details coincides with an even bigger enhancement in smooth regions. Indeed, with EOS-MD, ALVS and ALVD are increased w.r.t the scores of the input image by respectively ×2.73 and ×2.51. The proposed approach has lower ALVD compared to the other two methods. However, the smooth regions are much less enhanced too. Indeed, ALVS and ALVD are increased w.r.t the scores of the input image by respectively ×2.03 and ×2.18. Thus, thanks to the proposed algorithm we are able to enhance more the details than the noise. Finally, we did not notice halo artifacts on any outcomes of the three compared algorithms. However, the lowest ALVE value was obtained with the proposed method, which suggests that it may be less prone to this issue too.

CONCLUSIONS

In this paper we proposed an original method to decompose and enhance X-ray MD images that exploits the information provided by our adaptive NL means filter directly in the MS decomposition loop. An optimal representation of the details levels is obtained to limit the enhancement to relevant information of the image while avoiding too regular outcomes. The efficiency of the algorithm is demonstrated by quantitative and qualitative assessments on both phantoms and clini-cal MD images. As shown by the scores obtained on the phantoms, it seems possible to extend the method to LD images too. Thus, we intent to fulfill our study with tests on clinical LD images. Besides, we want to quantify how much we can decrease the amount of dose and still be able, thanks to the proposed method, to obtain a diagnostic image that targets other applications than follow-up exams. In this perspective, we are going to involve radiologists in the study to get a quantitative clinical assessment of the image quality.

Fig. 1 :

 1 Fig. 1: Architecture of the proposed method. I = input image; Î = denoised image; D i I = detail image at level i from the decomposition of the image I.

Fig. 2 :

 2 Fig. 2: Crop of the phantom 70kV and 200mA: (a) EOS-nD, (b) Proposed

Fig. 3 :

 3 Fig. 3: Crop from a spine image; dose = 14.64 µGy: (a) Low contrast and noisy image, (b) EOS-LD, (c) EOS-MD, (d) Our method.

  and 2. In details, with

	Params. kV mA C Dose µGy SNR CNR DYN SNR CNR DYN EOS-LD Prop.
	70 200 4	79.36	18.81 0.92 55% 35.05 3.70 61%
	90 200 4	124.04	26.92 1.57 51% 49.73 5.82 60%
	120 200 4	197.94	39.82 2.12 53% 67.21 7.30 61%

Table 1 :

 1 Quantitative assessment on phantoms: Low Dose

	Params. kV mA C Dose µGy SNR CNR DYN SNR CNR DYN EOS-MD Prop.
	60 80 3	10.29	12.93 2.2	62% 16.22 3.51 70%
	65 80 4	16.97	15.75 1.37 56% 22.31 2.58 63%
	72 125 4	34.08	19.29 2.30 58% 26.99 3.80 64%

Table 2 :

 2 Quantitative assessment on phantoms: Micro Dose2

Table 3 :

 3 Average ALV scores on the clinical data set

		ALVS	ALVD	ALVE
	Input	442.98	561.93	848.83
	EOS-LD	1175.18	1375.49	1792.57
	EOS-MD 1210.09	1411.26	1756.96
	Prop.	899.55	1227.70	1657.70

The image in the linear domain is obtained after the calibration that is the first step of the image processing reconstruction chain. Then, considering that the X-rays follow an exponential law, the image is converted to the logarithmic domain and we obtain I.