Mikael Carlavan
email: carlavan@i3s.unice.fr

Laure Blanc-Féraud

Marc Antonini

Carole Thiebaut
email: carole.thiebaut@cnes.fr

Christophe Latry
email: christophe.latry@cnes.fr

Yves Bobichon
email: yves.bobichon@thalesaleniaspace.com

Joint coding-denoising optimization of noisy images

Keywords: Global optimization, noisy source coding, quantizing, image denoising

In this paper, we propose to study the problem of noisy source coding/denoising. The challenge of this problem is that a global optimization is usually difficult to perform as the global fidelity criterion needs to be optimized in the same time over the sets of both coding and denoising parameters. Most of the bibliography in this domain is based on the fact that, for a specific criterion, the global optimization problem can be simply separated into two independent optimization problems: The noisy image should be first optimally denoised and this denoised image should then be optimally coded. In many applications however, the layout of the acquisition imaging chain is fixed and cannot be changed, that is a denoising step cannot be inserted before coding. For this reason, we are concerned here with the problem of global joint optimization in the case the denoising step is performed, as usual, after coding/decoding. In this configuration, we show how to express the global distortion as a function of the coding and denoising parameters. We present then an algorithm to minimize this distortion and to get the optimal values of these parameters. We show results of this joint optimization algorithm on classical test images and on a high dynamic range image, visually and in a rate-distortion sense.

I. INTRODUCTION

Images acquired by imaging systems are most of the time degraded by noise, which mainly comes from the imperfections of optical instruments. It is well-known that noise decreases the performances of coding schemes as it reduces the correlation between pixels [START_REF] Lo | Noise impact on error-free image compression[END_REF]. This problem is commonly referred as the noisy source coding problem [START_REF] Dobrushin | Tsybakov Information transmission with additional noise[END_REF]. Many works have been devoted to address this issue [START_REF] Ayanoglu | On optimal quantization of noisy sources[END_REF], [START_REF] Ephraim | A unified approach for encoding clean and noisy sources by means of waveform and autoregressive model vector quantization[END_REF], [START_REF] Sakrison | Source encoding in the presence of random disturbance[END_REF], [START_REF] Wolf | Transmission of noisy information to a noisy receiver with minimum distortion[END_REF], [START_REF] Rebollo-Monedero | Wyner-Ziv quantization and transform coding of noisy sources at high rates[END_REF], [START_REF] Fischer | Estimation and noisy source coding[END_REF] and [START_REF] Al-Shaykh | Lossy compression of noisy images[END_REF]. The majority of the mentioned works are based on the study of the global distortion optimization initially proposed by [START_REF] Sakrison | Source encoding in the presence of random disturbance[END_REF] and refined in [START_REF] Wolf | Transmission of noisy information to a noisy receiver with minimum distortion[END_REF]. This study states that the global distortion, if measured by the Mean Square Error (MSE), can be treated as two separated problems. First, the original source image should be optimally, in the minimum MSE sense, estimated from the noisy data and this estimate should then be optimally coded [START_REF] Wolf | Transmission of noisy information to a noisy receiver with minimum distortion[END_REF].

However, adding a supplementary step before coding is not always possible and in many cases the noisy acquired images are directly encoded without pre-processing. Critical applications such as satellite imaging can not indeed afford to insert pre-processing steps in the acquisition imaging chain as the onboard resources are highly limited. But one is still interested in optimizing the global imaging chain to obtain the best final image. So the imaging chain has to be considered as it is and the global distortion of this chain needs to be optimized. This is the focus of this paper.

More precisely, we are considering here the problem of optimal joint coding/denoising of a noisy image and focus on the acquisition imaging chain depicted figure 1. We show that, under certain hypotheses that we will describe, a closed-form expression of the global distortion can be obtained. We propose then to optimize this distortion, with respect to the coding and the denoising parameters, to reach the minimum global distortion. The originality of the proposed approach relies on the fact that we propose a global joint optimization which takes into account all the parameters of the imaging chain. We will also show that treating the optimization of coding/denoising as two separate problems (as in [START_REF] Rebollo-Monedero | Wyner-Ziv quantization and transform coding of noisy sources at high rates[END_REF], [START_REF] Fischer | Estimation and noisy source coding[END_REF], [START_REF] Al-Shaykh | Lossy compression of noisy images[END_REF]) is suboptimal when the denoising is performed after coding. More precisely, we emphasize the necessity to take into account the denoising step in the rate-distortion allocation of the coder. And, as we will see, this requirement is confirmed by the results which display a significant improvement in comparison to the classical method which executes coding and denoising separately.

The paper is organized as follows. In section II, we present the studied imaging chain and introduce January 10, 2013 DRAFT hypotheses and notations. We detail in section III the proposed approach and we show how to get a closed-form expression of the global distortion for the studied case. We detail the joint optimization of this distortion in section IV and we present the algorithm to get the optimal parameters of the coding and the denoising steps. We show results, visually and in a rate-distortion sense, of the proposed joint optimization algorithm on classical test images and on a remote sensing image, in section V. We conclude in section VI and present perspectives for future works.

II. HYPOTHESES AND NOTATIONS

In the following, the operators applied to the image are denoted with a bold uppercase letter. The non-bold uppercase letters represent random variables whose realizations are denoted by a lowercase letter. With this notation, x is a realization of the random variable X. (X) i denotes the ith element of the random variable X. These variables are multidimensional x ∈ R N where N is the number of pixels. W x is a random variable associated to the wavelet transform of x and we denote W x,j , j ∈ {0, . . . , J -1} (J being the number of subbands) the jth subband of the random variable W x . We have w x,j ∈ R Nj where N j is the size of the subband. Finally, we suppose that a wavelet subband w x,j follows a generalized centered Gaussian distribution law of parameter α wx,j > 0 and variance σ 2 wx,j > 0 [START_REF] Antonini | Image coding using wavelet transform[END_REF]. A wavelet subband probability density function p wx,j can then be modeled as

p wx,j (w x,j ) = A α wx,j σ wx,j e -B(αw x,j ) w x,j σw x,j αw x,j , (1)
with

A α wx,j = α wx,j B α wx,j 2Γ 1/α wx,j (2)
B α wx,j = Γ 3/α wx,j Γ 1/α wx,j , (3)
and Γ is the usual Gamma function. The parameters σ 2 wx,j and α wx,j of the distribution law will be estimated using the kurtosis-based technique proposed in [START_REF] Kasner | Universal trellis coded quantization[END_REF]. Note that the same assumption will be applied to all wavelet transforms in the chain with, of course, different distribution parameters. January 10, 2013 DRAFT

As mentioned previously, we study the imaging chain shown figure 1. We consider the special case of coding techniques based on wavelet transforms [START_REF] Shapiro | Embedded image coding using zerotrees of wavelet coefficients[END_REF], [START_REF] Said | A new, fast, and efficient image codec based on set partitioning in hierarchical trees[END_REF] and [START_REF] Taubman | High performance scalable image compression with EBCOT[END_REF]. The coding step is then approximately decomposed in a non-redundant wavelet transform followed by a scalar subband quantizer. Note that this approximation is actually close to the coding schemes presented in the cited works.

We also consider that the denoising step is performed in the same wavelet basis than the coding. This choice may however need further explanations. Usually, an efficient wavelet transform for image denoising strongly differs from a wavelet transform suited for image coding. Image denoising techniques actually require redundant wavelet transforms to represent the characteristics of an image such as contours and oriented details while increasing the number of coefficients in image compression may be problematic [START_REF] Chappelier | Oriented Wavelet Transform for Image Compression and Denoising[END_REF]. Hence, a non-redundant wavelet transform used for image compression leads most of the time to poor denoising results. We are however very confident that using the same basis for both coding and denoising may provide an optimized decoding-denoising structure gathered in a single fast and lowressources algorithm. Extending the current work to complex denoising schemes such as [START_REF] Donoho | De-noising by soft-thresholding[END_REF] is a difficult task that will be addressed in future works.

Based on these considerations, the studied imaging chain is represented in detail figure 2. In this chain, Fig. 2. Considered imaging chain we consider the instrumental noise z to be independent, identically distributed and to follow a centered normal distribution with variance σ 2 z . W is a wavelet transform, W its inverse. We denote w x,j and w z,j to be respectively the subband j of the wavelet transform of x and z. Each quantized subband w ỹ,j will be coded using an entropy encoder. As this operation does not introduce any degradation, it does not appear on the chain displayed figure 2. R is a linear restoration algorithm which operates independently on the wavelet coefficients of each subband j of the image and writes

w x,j = arg min w -w ỹ,j 2 2 + λ j w 2 2 subject to w ∈ R Nj , (4)
where λ j > 0 is a regularizing parameter. The restoration algorithm (4) has a closed-form solution which January 10, 2013 DRAFT writes

w x,j = w ỹ,j 1 + λ j . (5)
We are aware of the simplicity of the considered algorithm, it appears however that the linearity of the restoration algorithm R is required if one wants to write the global distortion in closed-form. As mentioned previously, much work need to be addressed to consider state-of-art denoising algorithms.

The quantizer Q is an infinite mid-tread scalar subband quantizer of step ∆ j > 0 and is modeled as

Q(w y,j ) = ∆ j w y,j ∆ j + 1 2 , (6)
where ⌊ ⌋ is the floor function which returns the greatest integer less than or equal to its argument. We now detail the basis of the proposed method to compute a closed-form expression of the global distortion.

Let w b,j be the coding error of the subband j w b,j = Q(w y,j ) -w y,j .

We have

w ỹ,j = Q(w y,j ) = w y,j + w b,j = w x,j + w z,j + w b,j = w x,j + w ǫ,j , (8)
where w ǫ,j = w z,j + w b,j . The main hypothesis of the proposed method is to consider the first-order moment of the term w ǫ,j to be independent to the one of w x,j , that is

E [W ǫ,j W x,j ] = E [W ǫ,j ] E [W x,j ] , (9)
where W ǫ,j and W x,j are the random variables associated to w ǫ,j and w x,j . This hypothesis is mainly based on the fact that the quantizing part of the scheme figure 2 can be seen as a non-substractive dithering system where the Gaussian instrumental noise acts as a dithering noise.

More generally, a dithering system consists in inserting a noise with a certain probability density function prior to quantizing, to improve the decorrelation property [START_REF] Vanderkooy | Dither in Digital Audio[END_REF]. As mentioned in [START_REF] Wannamaker | A theory of nonsubtractive dither[END_REF], a non-substractive dithering system (named non-substractive as the dithering noise is not substracted after quantizing) allows the moments of the global error (that is the sum of the coding error and dithering noise) to be fully decorrelated to the moments of the coding source.

It happens that a Gaussian distribution stands among the probability density functions which allow a January 10, 2013 DRAFT noise to be considered as a dithering noise. The idea here is then to take benefit of the presence of the instrumental noise by considering it as a dithering noise. With such consideration, we know that the m first-order moments of the global error are deccorelated to the n first-order moments of the quantizing source. In detail, we have for any integer m > 0 and n > 0

[17] E W m ǫ,j W n x,j = E W m ǫ,j E W n x,j (10)
We focus on the special case m = 1 and n = 1, giving the propery ( 9), as it is the basis of the proposed approach. Moreover, if the instrumental noise z meets the dithering noise requirements, we also have

[17] E [W ǫ,j ] = 0, (11)
E W ǫ,j 2 = N j σ 2 wz,j + N j ∆ 2 j 12 , (12)
where σ wz,j is the standard deviation of the distribution law of the wavelet transform w z,j . From [START_REF] Vanderkooy | Dither in Digital Audio[END_REF] we know that a Gaussian noise effectively owns the properties of a dither noise if the standard deviation of its distribution law is large enough. In the present case, the condition (10) will be verified if the following statement is true

σ wz,j > ∆ j 2 . ( 13
)
As the standard deviation of instrumental noise is usually low in imaging systems, the condition [START_REF] Said | A new, fast, and efficient image codec based on set partitioning in hierarchical trees[END_REF] assumes that the proposed approach will be valid only for asymptotic coding rate. We will however develop our method to consider all coding rates.

III. GLOBAL RATE-DISTORTION ANALYSIS

As mentioned in the section II, the studied imaging chain depends on two sets of parameters: The regularizing parameters λ j in (5) and the quantizing steps ∆ j in ( 6), for each j ∈ {0, . . . , J -1}. The global rate-distortion joint optimization problem consists in finding the optimal parameters λ * j and ∆ * j which minimize the global distortion D under the constraint that the coding rate R does not exceed the target rate R c . This can be formalized as the following

λ * j , ∆ * j = arg min D(λ j , ∆ j ) subject to R(λ j , ∆ j ) ≤ R c , λ j > 0, ∆ j > 0 . ( 14
)
Under this form, the optimization problem ( 14) is difficult to solve so that it is usually written under an unconstrained form. Let τ > 0 be a Lagrange mutliplier. Problem [START_REF] Taubman | High performance scalable image compression with EBCOT[END_REF] can then be written [START_REF] Everett | Generalized Lagrange Multiplier Method for Solving Problems of Optimum Allocation of Resources[END_REF] λ

* j , ∆ * j = arg min D(λ j , ∆ j ) + τ (R(λ j , ∆ j ) -R c ) subject to λ j > 0, ∆ j > 0, τ > 0 . ( 15
)
To solve the global distortion joint optimization problem [START_REF] Chappelier | Oriented Wavelet Transform for Image Compression and Denoising[END_REF], we then need to express the global distortion D and the global coding rate R as a function of the regularizing parameters λ j and the quantizing steps

∆ j .
Proposition 1: If σ wz,j verifies hypothesis [START_REF] Said | A new, fast, and efficient image codec based on set partitioning in hierarchical trees[END_REF] for each j ∈ {0, . . . , J -1}, then the global distortion D of the imaging chain displayed figure 2 writes

D = J-1 j=0 π j a j λ 2 j (1 + λ j ) 2 σ 2 wx,j + π j a j (1 + λ j ) 2 σ 2 wz,j + π j a j (1 + λ j ) 2 ∆ 2 j 12 , (16)
where

a j = N j N , (17)
is the weight of the subband j in the whole image.

Proof: We start from the fact that the (mean) global distortion writes

D = 1 N E X -X 2 , ( 18
)
where X is the random variable associated to the output final image x. Thanks to the orthogonality of the wavelet subbands, the global distortion can also be formulated as

D = 1 N J-1 j=0 π j E W x,j -W x,j 2 , (19)
where π j are weighting coefficients which depend on the filters and the decimation factors used in the wavelet transform [START_REF] Usevitch | Optimal bit allocation for biorthogonal wavelet coding[END_REF]. Note that these weighting coefficients are only required if one considers biorthogonal wavelet transforms such as the Cohen-Daubechies-Feauveau (CDF) 9/7 wavelet transform [START_REF] Cohen | Biorthogonal bases of compactly supported wavelets[END_REF]. They are equal to 1 for an orthogonal wavelet transform.

In the case of the studied imaging chain displayed figure 2, the final image is the output of the restoration and writes

w x,j = Rw ỹ,j . (20)
Using equations ( 5) and ( 8), the final image can be expressed as a function of the source and the global error

w x,j = w x,j 1 + λ j + w ǫ,j 1 + λ j . (21)
From equations ( 19), ( 21) and using the moments decorrelation hypothesis (9), we deduce the global distortion

D = 1 N E X -X 2 = 1 N J-1 j=0 π j λ 2 j (1 + λ j ) 2 E W x,j 2 + π j (1 + λ j ) 2 E W ǫ,j 2 . ( 22
)
Finally, the global distortion ( 22) can be further developed using the results [START_REF] Shapiro | Embedded image coding using zerotrees of wavelet coefficients[END_REF] to obtain the expression [START_REF] Donoho | De-noising by soft-thresholding[END_REF].

Note that the global distortion ( 16) requires the knowledge of the variance of each subband of the original image σ 2 wx,j . This variance is generally unknown but can be deduced from the observed image if we consider an orthogonal wavelet transform. In that case, the variance of the noise in each wavelet subband j is equal to the variance of the noise in the image domain, i.e. σ 2 wz,j = σ 2 z , which is supposed to be known. Then, σ 2 wx,j can be computed during the rate-allocation of the coder from the observed subband variance σ 2 wy,j using the fact that

σ 2 wx,j = σ 2 wy,j -σ 2 z . (23)
The second part of the problem [START_REF] Chappelier | Oriented Wavelet Transform for Image Compression and Denoising[END_REF] requires the expression of the global coding rate R. This rate can be expressed as the weighted sum of the rate in each subband

R j R = J-1 j=0 a j R j (∆ j ), (24)
where a j is given in [START_REF] Wannamaker | A theory of nonsubtractive dither[END_REF]. As mentioned in the hypotheses section, we assume that each quantized subband is encoded using an entropy encoder. Then, in a general non-asymptotic case, the coding rate of a subband j can be estimated by its entropy [START_REF] Shannon | A Mathematical Theory of Communication[END_REF] R j (∆ j ) = -+∞ m=-∞ P wy,j (m, ∆ j ) log 2 P wy,j (m, ∆ j ) ,

where P wy,j (m, ∆ j ) is the probability to get the symbol m which depends on the density probability function p wy,j , defined in (1), of the subband w y,j and on the quantizing step ∆ j P wy,j (m, ∆ j ) =

m∆j+ ∆ j 2 m∆j- ∆ j 2
p wy,j (w y,j )dw y,j .

January 10, 2013 DRAFT

As mentioned in the hypothesis section, we assume that each wavelet subband followed the generalized centered Gaussian distribution law defined in [START_REF] Lo | Noise impact on error-free image compression[END_REF]. The density probability function p wy,j is then given by: p wy,j (w y,j ) = A α wy,j σ wy,j e -B(αw y,j )

w y,j σw y,j αw y,j ,

A α wy,j = α wy,j B α wy,j 2Γ 1/α wy,j

B α wy,j = Γ 3/α wy,j Γ 1/α wy,j ,

and where σ 2 wx,j and α wx,j are the parameters of the distribution law, estimated using the kurtosis-based technique detailled in [START_REF] Kasner | Universal trellis coded quantization[END_REF].

Proposition 2: The global rate-distortion optimization problem ( 14) can be solved by minimizing

φ(∆ j , λ j , τ ) = J-1 j=0 π j a j λ 2 j (1 + λ j ) 2 σ 2 wx,j + π j a j (1 + λ j ) 2 σ 2 wz,j + π j a j ∆ 2 j 12(1 + λ j ) 2 + τ   J-1 j=0 a j R j (∆ j ) -R c   , (30)
with respect to ∆ j > 0, λ j > 0 and τ > 0; R j (∆ j ) being defined in [START_REF] Shannon | Coding theorems for a discrete source with a fidelity criterion[END_REF].

Proof: This demonstration is straightforward. From equation ( 15), we define

φ(∆ j , λ j , τ ) = D + τ (R -R c ) , (31)
and we substitute D and R with their respective expressions ( 16) and [START_REF] Parisot | Optimization of the Joint Coding/decoding Structure[END_REF].

We detail in the next part how to minimize the criterion (30).

IV. GLOBAL RATE-DISTORTION OPTIMIZATION

Using proposition 2, the optimization problem (15) becomes

∆ * j , λ * j = arg min φ(∆ j , λ j , τ ) subject to λ j > 0, ∆ j > 0, τ > 0 . ( 32
)
The convexity analysis of the function φ is a difficult task, such that the existence and uniqueness of solutions of problem (32) is not straightforward. We explicitly assume that a minimum of problem (32) exists, and we will show that among all the stationary points of the function φ, only one meets the positivity constraint of the parameters. We propose a numerical algorithm to find this minimum. This algorithm is based on the resolution of the simultaneous equations obtained from the Karush-Kuhn-Tucker January 10, 2013 DRAFT (KKT) conditions [START_REF] Kuhn | Nonlinear programming[END_REF] of problem (32).

Proposition 3: The KKT conditions of problem (32) admits only one solution (λ * j , τ * , ∆ * j ) which verifies

λ * j = σ 2 wz,j σ 2 wx,j + ∆ * j 2 12σ 2 wx,j . ( 33
)
π j ∆ * j 6(1 + λ j ) 2 + τ * ∂R j ∂∆ j (∆ * j ) = 0. ( 34
) J-1 j=0 a j R j (∆ * j ) = R c (35)
Proof: From the KKT conditions of problem (32), we get

               ∂φ(∆ * j ,λ * j ,τ * ) ∂∆j = ajπj∆ * j 6(1+λ * j ) 2 + τ * a j ∂Rj ∂∆j (∆ * j ) = 0 ∂φ(∆ * j ,λ * j ,τ * ) ∂τ = J-1 j=0 a j R j (∆ * j ) -R c = 0 ∂φ(∆ * j ,λ * j ,τ * ) ∂λj = 12ajπjλ * j σ 2 w x,j -12ajπjσ 2 w z,j -ajπj∆ * j 2 6(1+λ * j ) 3 = 0 (36)
with

∂R j ∂∆ j (∆ j ) = - 1 log(2) +∞ m=-∞ 1 + log P wy,j (m, ∆ j ) × p wy,j m∆ j + ∆ j 2 m + 1 2 -p wy,j m∆ j - ∆ j 2 m - 1 2 . ( 37
)
The expressions (33), (34) and (35) of the optimal parameters directly follow from the optimality conditions (36). The uniqueness of τ * can be proved as follows. Using equations ( 33) and (34), for a given τ > 0, we define

g τ (∆) = π j ∆ j 6 1 + σ 2 w z,j σ 2 w x,j + ∆j 2 12σ 2 w x,j 2 + τ ∂R j ∂∆ j (∆ j ). (38)
The coding rate R j is a monotonically decreasing positive function with respect to ∆ j [START_REF] Shannon | Coding theorems for a discrete source with a fidelity criterion[END_REF], ∆ j being positive. Its limits are zero when ∆ j tends to infinity and infinity when ∆ j vanishes to zero [START_REF] Gish | Asymptotically efficient quantizing[END_REF]. Its derivative ∂Rj ∂∆j is negative and monotonically increasing, whose limits are minus infinity when ∆ j vanishes to zero and zero when ∆ j tends to infinity [START_REF] Shannon | Coding theorems for a discrete source with a fidelity criterion[END_REF].

We can see that the function g τ tends to minus infinity when ∆ j vanishes to zero. The limit of g τ when ∆ j tends to infinity is more difficult to address. When ∆ j tends to infinity, the first term of g τ is positive and evolves as 1 ∆ 3 j . The second term is negative but depending on the value of τ , the function g τ January 10, 2013 DRAFT may not cross zero such that the optimal quantizing step ∆ * j (which verifies g τ (∆ * ) = 0) may not always exist. However, this never happened in our simulations and we found out that the term ∂Rj ∂∆j (∆ j ) tends to zero quickly, when ∆ j tends to infinity, such that g τ always crosses zero and therefore ∆ * j always exists. The uniqueness of ∆ * j depends on the number of times the function g τ crosses zero. Using the first derivative of the function g τ , we can show that g τ crosses zero only once. We have

g ′ τ (∆) = π j 6 1 + σ 2 w z,j σ 2 w x,j + ∆j 2 12σ 2 w x,j 3 1 + σ 2 wz,j σ 2 wx,j - ∆ 2 j 4σ 2 wx,j + τ ∂ 2 R j ∂∆ 2 j (∆ j ), (39)
To prove that g τ crosses zero only once, we only study the sign of the first derivative (39). From [START_REF] Shannon | Coding theorems for a discrete source with a fidelity criterion[END_REF],

we can say that ∂ 2 Rj ∂∆ 2 j is positive and monotonically decreasing. From (39), we see that g ′ τ is positive when ∆ j vanishes to zero. When ∆ j tends to infinity, the first term of g ′ τ evolves as -1 ∆ 4 . We use the same remark that the one we used for the study of g τ limits and we assume that ∂ 2 Rj ∂∆ 2 j decreases very quickly, when ∆ j tends to infinity, such that g ′ τ reduces to its first term. In that case, we can say that g ′ τ is negative when ∆ j tends to infinity and, therefore, changes its sign only once.

From this result, we can conclude that the function g τ starts from minus infinity, crosses zero to reach its maximum and then decreases back to zero. This function crosses zero only once, implying that the optimal quantizing step ∆ * j , if it exists, is unique for any given τ > 0. In this paper, we propose to compute the value of ∆ * j numerically using a binary search procedure. But many precautions need to be taken to be sure that we always operate on the increasing part of the function g τ .

From (38), we see that ∆ * j is function of τ . It seems reasonable to think that the higher τ is, the higher ∆ * j needs to be for the function (38) to cross zero. This implies that the optimal quantizing step ∆ * j can then be noted as a function of τ

∆ * j = f (τ ), (40)
where f is an increasing function. Consequently, from [START_REF] Shannon | Coding theorems for a discrete source with a fidelity criterion[END_REF], we deduce that the coding rate R j is a monotonically decreasing function with respect to τ . Using (35) and (40), we define

h(τ ) = J-1 j=0 a j R j (f (τ )) -R c . (41)
Then it seems clear that the function h is a monotonically decreasing function with respect to τ whose limits are infinity when τ vanishes to zero and -R c when τ tends to infinity. Its root τ * , which verifies h(τ * ) = 0, exists, is unique and can be computed using any root-finding algorithm. In our simulations, a binary search procedure will be used.

We previously showed that for any given τ > 0, the optimal quantizing step ∆ * j exists and is unique. Using the uniqueness of τ * we can conclude that the parameters τ * and ∆ * j which verify (35) and (34) are unique. From this result and equation (33), the existence and uniqueness of the optimal regularizing parameter λ * j is straightforward to show. For the sake of simplicity, each binary search algorithm will be parametrized to the same given precision ρ = 0.1. The case of the low frequency subband (j = J -1) will be processed differently as we do not want to degrade these coefficients. We will only use quantizing to round these coefficients to their nearest integers. Consequently, we will set

∆ * J-1 = 1, (42)
λ * J-1 = σ 2 wz,J-1 σ 2 wx,J-1 + 1 12σ 2 wx,J-1 . (43)
Finally, the overall joint optimization procedure for solving problem ( 14) is given in the algorithm 1.

Note that the binary search sub-procedures are not detailled in this process. The algorithm 1 intends to be quite general and we let the choice of the root-finding algorithms to the user.

V. RESULTS

We simulate the joint optimization algorithm 1 on the well-known test images Lena, Barbara, Pirate and on the high-dynamic range remote sensing image displayed figure 3. For this simulation, we set the wavelet transform W to be a three levels CDF 9/7 wavelet transform [START_REF] Cohen | Biorthogonal bases of compactly supported wavelets[END_REF] and R is given by (4). Each test image has been noised with an additive white Gaussian noise with a standard deviation σ z equal to 15. As the efficiency of the proposed estimation depends on the standard deviation of this noise (see equation ( 13)), we simulate two more cases for the Barbara image: σ z = 5 and σ z = 25.

For each target rate, we simulate the imaging chain 2 with the usual disjoint optimization technique, which consists in selecting the quantizing steps and the regularizing parameters such that the coding and the restoration errors are independently minimized. The coding error minimization has been achieved using the rate-distortion allocation based model proposed in [START_REF] Parisot | Optimization of the Joint Coding/decoding Structure[END_REF]. As for the restoration error, it has been minimized using an exhaustive search of the optimal regularizing parameters. Once the final image has been reconstructed using these parameters, we numerically compute the global distortion

D = 1 N x -x 2 , ( 44
)
where x is the clean test image (assumed to be known in our numerical experiments) and x is the final image. The distortion (44) is the true distortion and will be referred as the ground truth in our simulations. fact that for low coding rates, the condition ( 13) is not respected anymore and that the moments of the global error cannot consequently be considered decorrelated to the moments of the source. As mentioned previously, we performed several simulations on the Barbara image with different standard deviations of the noise such that the condition (13) can be verified for different ranges of coding rate. When the standard deviation is low (figure 5), we see that the proposed estimation is performant if the coding rate is around 2.5 bits/pixel and more. However for this high coding rate, the coding step is almost lossless such that the global optimization problem is reduced to the optimization of the restoration only.

Therefore, the joint and the disjoint optimization techniques become the same and give then similar results.

January 10, 2013 DRAFT As shown by the figures 5, 6 and 7, the range of validity of the proposed estimation increases as the standard deviation increases. For a high standard deviation (figure 7), we can verify that the proposed estimation is valid for lower coding rates (around 1.8 bits/pixel and more). In that case, the joint optimization displays significant improvement in comparison to the disjoint optimization. It allows for example to reach the same global error than the disjoint optimized technique but for a lower coding rate. On the Barbara image and for σ z = 15 (figure 6), the joint optimization technique reaches at 1.42 bits/pixel the same distortion than the one obtained at 2.04 bits/pixels for the disjoint optimization technique, saving therefore almost 30% of the bit budget. The benefit in term of compression performances of the joint optimization technique appears then to be very significant.

Visual results for the target rate of 2.5 bits/pixel are given figures 9 to 12 for the high-dynamic range remote sensing image. We do not focus on the quality of the reconstructed images regarding to the reference one as the considered chain is excessively simple. Clearly, the presence of artifacts on the reconstructed image is due to the simple hypothesis that we made on the restoration algorithm (4). On the contrary, we are more concerned on the improvement of the image quality of the joint optimized chain with respect to the disjoint optimized one. We can see that the global joint optimization of the chain always leads to a reconstructed image witch contains less blurry edges or ringing artifacts. This is particularly visible on the edges of the buildings figures 9 and 11. As mentioned in the introduction of this paper, the obtained results clearly point that optimizing coding and denoising separately is suboptimal.

One needs instead to address the problem of imaging chain design in its globality; the proposed method and the obtained results are encouraging in this sense. A lot of works is however required to extend the proposed method to lower coding rates and to more complex denoising schemes.

VI. CONCLUSIONS

In this paper we considered the problem of joint noisy source coding/denoising. Most of the time, the coding and the denoising parameters are selected independently such that the coding and the restoration error are respectively minimized. This parameters selection technique leads however to a suboptimal distortion. It appears then crucial to address the problem of joint coding/denoising in its globality. We proposed here a technique to modelize the global distortion and we presented an algorithm to get the optimal coding and denoising parameters. We simulated this joint optimization technique on classical test January 10, 2013 DRAFT

Fig. 1 .

 1 Fig. 1. Considered acquisition imaging chain.

Fig. 3 .

 3 Fig. 3. Test images: (a) is Lena (256 × 256 pixels), (b) is Barbara (512 × 512 pixels), (c) is Pirate (1024 × 1024 pixels) and (d) is a high-dynamic range remote sensing image (12 bits) of Cannes harbour (1024 × 1024 pixels).

Fig. 4 .Fig. 5 .

 45 Fig. 4. Comparison of the disjoint optimized distortion (ground truth and model-based estimation) to the joint optimized distortion (model-based estimation) on Lena, σz = 15.

Fig. 6 .Fig. 7 .

 67 Fig. 6. Comparison of the disjoint optimized distortion (ground truth and model-based estimation) to the joint optimized distortion (model-based estimation) on Barbara, σz = 15.

Fig. 8 .

 8 Fig. 8. Comparison of the disjoint optimized distortion (ground truth and model-based estimation) to the joint optimized distortion (model-based estimation) on Pirate, σz = 15.

Fig. 9 .Fig. 10 .

 910 Fig. 9. Visual comparison of reconstruction results. (a) is the reference image, (b) is the observed image, (c) is the image reconstructed with the parameters obtained by the minimizing the measured distortion and (d) is the image reconstructed with the parameters obtained by the model-based optimization algorithm 1. The coding rate is 2.5 bits/pixel. The image range has been extended to point up the image reconstruction artifacts.

Fig. 11 .Fig. 12 .

 1112 Fig. 11. Visual comparison of reconstruction results. (a) is the reference image, (b) is the observed image, (c) is the image reconstructed with the parameters obtained by the minimizing the measured distortion and (d) is the image reconstructed with the parameters obtained by the model-based optimization algorithm 1. The coding rate is 2.5 bits/pixel. The image range has been extended to point up the image reconstruction artifacts.

January 10, 2013 DRAFT

images and on a high-dynamic range remote sensing image. We concluded that our joint coding/denoising optimization approach can either allows to reach the same quality at lower rates or to improve the quality of the reconstructed final image for the same rates, in comparison to the image obtained using the classical disjoint optimization technique. Further works will be focussed on the extension of the proposed model to lower coding rates and to advanced denoising schemes.

Algorithm 1 Global rate-distortion joint optimization algorithm

Set τ = 1.

Set ρ = 0.1.

while

J-1 j=0 a j R j -R c > ρ do for j from 0 to J -2 do

Compute the value of the regularizing parameter λ j from (33).

Compute the value of the regularizing parameter λ j from (33).

end while end for

Compute the regularizing paramater λ J-1 from (42).

Compute the quantizing step ∆

end if end while

Output the optimal regularizing parameters λ * j . Output the optimal quantizing steps ∆ * j .

The estimation model [START_REF] Donoho | De-noising by soft-thresholding[END_REF] of the global distortion that we proposed has then been computed with the values of parameters obtained for the ground truth. This allows to verify that the estimation (16) of the global distortion is close to the ground truth (44), implying the validity of the proposed method. And finally, we use the proposed joint optimization algorithm 1 to compute the optimal parameters, that we inserted into the estimation model ( 16) to compute the minimal distortion.

The resulting rate-distortion curves are given figures 4 to 8. We see that the validity of the proposed estimation, as expected by the hypothesis [START_REF] Said | A new, fast, and efficient image codec based on set partitioning in hierarchical trees[END_REF], is not always verified and depends on the target coding rate. More precisely, the proposed estimation approximates well the true distortion for medium to high coding rates but does not give satisfying results for low coding rates. This can be explained by the January 10, 2013 DRAFT


Mikael Carlavan
email: carlavan@i3s.unice.fr

Laure Blanc-Féraud

Marc Antonini

Carole Thiebaut
email: carole.thiebaut@cnes.fr

Christophe Latry
email: christophe.latry@cnes.fr

Yves Bobichon
email: yves.bobichon@thalesaleniaspace.com

On the optimization of the satellite imaging chain

Keywords: Global optimization, satellite imaging chain, coding, image restoration I

In this paper, we focus on the global optimization of the satellite imaging chain. The theoretical analysis of the satellite imaging chain optimization is a difficult problem that needs lot of approximations.

In order to consider the complex real satellite imaging chain, we propose to address this problem numerically and we present, based on numerical experiments, techniques to optimize the quality of the reconstructed final image. We first focus on the common question of the position of the restoration step in the imaging chain, that is on-board before coding or on-ground after coding. Then, we present several methods to remove the coding artifacts inherent in wavelet based coder schemes. From these numerical results we propose a new satellite imaging chain and we show visual and rate-distortion results on a real satellite image.

performed on board. As shown on the figure 1, a satellite imaging chain is mainly composed by three parts: The acquisition which captures and samples the scene, the compression which reduces the volume of data of the image before storing it in the on-board memory (compression is also required for a fast transmission) and the restoration which attenuates the degradations collected during the processing chain.

For the rest of this paper, we assume that the transmission does not alter the quality of the image. Usually, the acquisition is fixed and fully characterized by the optics of the satellite and the electronics of the chain. The parameters of this part are therefore quite fixed such that the acquisition cannot be tuned so much to improve the quality of the final image. On the contrary, the compression and the restoration have many more degrees of freedom. The choice of the parameters of these two parts is therefore of crucial importance in the control of the satellite imaging chain. This is the problem of global optimization.

Formally, the global optimization of a satellite imaging chain is the problem of designing the optimal January 9, 2013 DRAFT coding and restoration algorithms which give the best final image knowing the parameters of the chain (sampling, blur, noise, coding rate, ...). This problem is very challenging and has been little investigated so far. Some works have been done in designing an optimal coding/decoding structure [START_REF] Parisot | Optimization of the Joint Coding/decoding Structure[END_REF] which takes into account the characteristics of the imaging chain, or an efficient restoration method which deals with the quantizing noise [START_REF] Tramini | Optimal joint decoding/deblurring method for optical images[END_REF]. But to the best of our knowledge, the study of the global system optimization, which includes both compression and restoration, has not devoted much work.

To formulate this specific problem, we consider the imaging chain showed figure 1

where R c is the target coding rate and E is the expectation operator with respect to the distribution law of x. Solving problem (1) is very difficult in many aspects. Firstly, problem (1) searches for the optimal coder and restoration among all available techniques, which is not tractable. Second, even if the coding and restoration methods are given and perfectly known, an analytic expression of the global distortion is usually not available and depends on the distance measure D and on the knowledge of the real unknown image x (or its statistics). As we can see, the problem (1) is difficult to solve in a general context. The main contribution of this paper is to show that this problem can however be simplified in the specific case of D being the well-known Mean Square Error (MSE).

From the results of [START_REF] Wolf | Transmission of noisy information to a noisy receiver with minimum distortion[END_REF], we show here that, in the case of the MSE criterion, the global joint-optimization problem (1) can be replaced by two disjoint optimization problems. More precisely, the original image

x should first be optimally estimated from the intrumental image y and, second, this estimation should be optimally coded. This result is particularly interesting in the case of satellite imaging as it addresses the position of the estimation process in the chain, that is the position of the restoration with respect to January 9, 2013 DRAFT the coding step. We will present in this paper numerical experiments on this aspect.

One difficulty of the optimization of the imaging chain is also the presence of coding artifacts in the reconstructed image. We also address this issue in the paper and we present techniques to treat this problem of coding noise removal. Finally, the last contribution of the paper is to present a new satellite imaging chain. The particularity of this chain is that it deduces from the conclusions obtained on each problematic (i.e. the position of the restoration and how to deal with the coding noise) and is therefore optimized on these aspects.

The paper is organized as follows. We present, in sections II and III, numerical experiments to improve the quality of the final image. The purpose of these parts is to bring leads, mainly based on numerical results, to open problems in the design of a satellite imaging chain such as the position of the restoration in the chain and how to deal with the coding noise. From the obtained results we propose a new imaging chain in section IV and we present reconstruction results on a real satellite image. Finally, section V concludes the paper and presents perspectives for future works.

The data presented in this paper are provided by the French Space Agence (CNES 1 ) and are simulations of the PLEIADES-HR satellite. We will then focus only on the imaging chain of this satellite but the methods we propose are more general and can be easily extended to the characteristics of other satellites.

II. THE POSITION OF THE RESTORATION IN THE IMAGING CHAIN

The actual restoration method used by the CNES only deals with the blur and the additive Gaussian noise of the instrument [START_REF] Lier | Imagerie spatiale : des principes d'acquisition au traitement des images optiques pour l'observation de la Terre[END_REF]. It actually does not take into account the fact that the decoded image is also deteriorated with coding noise. This restoration is therefore also suitable to be used on-board just before coding on the intrumental image, as this image perfectly matches the image formation model considered by the restoration.

In this part, we are thus considering the question of the position of the restoration in the satellite imaging chain, i.e. on-board 2 before coding (see figure 2) or as usual on-ground after coding (see figure 1).

1 Centre National d'Etudes Spatiales 2 The purpose of this part is to give a methodological answer to the position of the restoration in the chain; so we assume that we have at our disposal enough on-board computing resources such that an on-board restoration can be implemented

A. Introducing the on-board restoration

As mentioned in the introduction, the initial global optimization problem consists in finding the optimal coding C * and restoration T * such that

C * , T * = arg min E [D(x, T (C(y)))] subject to C, T R(C(y)) ≤ R c . (2)
A similar problem has been formalized in [START_REF] Wolf | Transmission of noisy information to a noisy receiver with minimum distortion[END_REF] for the study of optimal noisy source coding. The main result of [START_REF] Wolf | Transmission of noisy information to a noisy receiver with minimum distortion[END_REF] states that the global distortion, if measured by the Mean Square Error (MSE), can be separated in two terms. More precisely, the authors of [START_REF] Wolf | Transmission of noisy information to a noisy receiver with minimum distortion[END_REF] showed that

E x -T (C(y)) 2 2 = E x -E[x|y] 2 2 + E E[x|y] -T (C(y)) 2 2 , (3)
where

D * = E x -E[x|y] 2 2 + min E E[x|y] -T (C(y)) 2 2 subject to C, T . (4)
From ( 4), we propose to use the MSE as the global measure criterion and we rewrite the optimization problem (2) as

C * , T * = arg min E E[x|y] -T (C(y)) 2 2 subject to C, T R(C(y)) ≤ R c . ( 5
)
We see that the global distortion can be expressed and optimized with respect to the image E[x|y] instead of the original image x. Note that the problem ( 5) is not simpler to solve as the computation of the image E[x|y] is usually not accessible. This formulation may however be interesting to address the common question of the position of the restoration [START_REF] Lambert-Nebout | On-board optical image compression for future high-resolution remote sensing systems[END_REF] in the design of the satellite imaging chain (that is before or after coding).

As mentioned previously, the image E[x|y] represents the restoration of the true image x from the instrumental one y. It is then very tempting to think that this ideal image is actually the result of the restoration T , moved on-board of the satellite (i.e. before coding), and that we can replace E[x|y] in [START_REF] Lambert-Nebout | On-board optical image compression for future high-resolution remote sensing systems[END_REF] by T (y) such that the global optimization problem can also be written

C * , T * = arg min E T (y) -C(T (y)) 2 2 subject to C, T R(C(T (y))) ≤ R c . (6)
It is clear that the problem ( 6) is slightly different from the initial optimization problem (2) but problem ( 6) is easier to treat as each variable can almost be optimized separately. If T is fixed, then problem [START_REF] Shannon | A Mathematical Theory of Communication[END_REF] looks for the optimal coder C * which minimizes the coding error under the constraint that the coding rate does not exceed the target coding rate. This problem is well-known and referred as the coding rate-allocation problem [START_REF] Shannon | A Mathematical Theory of Communication[END_REF] which has been addressed a lot in the coding community [START_REF] Antonini | Image coding using wavelet transform[END_REF], [START_REF] Berger | Rate Distortion Theory: A Mathematical Basis for Data Compression[END_REF], [START_REF] Ortega | Rate-distortion methods for image and video compression[END_REF] and references therein.

To be clear, the global joint optimization problem (2) is very difficult to address. But, in our opinion, we believe that moving the restoration on-board could take benefit of the characteristics of the actual restoration and coding algorithms. So one way (but again this is not the only one) to address the problem As mentioned previously, the latter has been the focus of intense work in the imaging community. In the next part, we simulate the comparison between the on-board chain and the on-ground chain with the state-of-the-art restoration algorithms.

B. Comparison of on-board and on-ground chains

We are considering the on-board chain showed figure 2 in comparison to the on-ground one illustrated figure 1. For the simulation, the coding step is very similar to the one implemented in the PLEIADES-HR satellite and the restoration will be performed in two steps as follows. First, a direct deconvolution will be performed using the target Point Spread Function (PSF) provided by the CNES. The second step is the denoising step. A wavelet packet decomposition [START_REF] Kalifa | Deconvolution by thresholding in mirror wavelet bases[END_REF] is usually employed to fit the frequential characteristics of the deconvolved noise [START_REF] Lier | Imagerie spatiale : des principes d'acquisition au traitement des images optiques pour l'observation de la Terre[END_REF]. However, another important point to take into account for an efficient denoising is the decrease rate of reconstruction error from the M largest wavelet coefficients [START_REF] Patel | Shearlet-Based Deconvolution[END_REF]. The faster the reconstruction error decreases, the better the denoising is. And on this point, a wavelet packet transform is not optimal [START_REF] Mallat | A wavelet tour of signal processing[END_REF].

We propose here to perform the denoising using a variant of the wavelet transform named the Shearlet transform [START_REF] Labate | Sparse multidimensional representation using shearlets Proc[END_REF]. This recent transform is an extension of the wavelet transform and considers a "shear" matrix instead of the classical dyadic shifts and dilations. This offers the ability to capture oriented details and is, among the contourlets [START_REF] Do | The contourlet transform: an efficient directional multiresolution image representation[END_REF] and the curvelets [START_REF] Candès | Fast discrete curvelet transforms[END_REF], an optimal transform (in term of coefficients magnitude decreasing rate) for the representation of images [START_REF] Patel | Shearlet-Based Deconvolution[END_REF]. We will compare this method to the current state-of-art restoration methods such as the FORWARD method [START_REF] Neelamani | ForWaRD: Fourier-Wavelet Regularized Deconvolution for Ill-Conditioned Systems[END_REF] or the method based on a Stein block thresholding [START_REF] Chesneau | Stein block thresholding for wavelet-based image deconvolution[END_REF]. Note that to be coherent with the current technique used for the PLEIADES-HR satellite, we only focus here on the methods which decompose the restoration in a direct deconvolution followed by a threshold operation of some sparse representation. We did not include the methods based on a variational framework such as [START_REF] Bect | A l1-unified variational framework for image restoration[END_REF].

We simulate these two chains on the image presented figure 3 and we compare these results to the reconstructed images provided by the CNES. From numerical experiments, we observed that the threshold parameters which minimize the MSE almost follow a power of two decreasing law (i.e. at each level, the thresholds parameters are divided by 2). The threshold parameters of each method have then been computed using this technique. In this simulation, the original image x is known and the MSE can thus be computed. Note that in a real environnment, unbiased estimator the MSE exist and do not require the knowledge of the true image [START_REF] Ramani | Monte-Carlo SURE: A Black-Box Optimization of Regularization Parameters for General Denoising Algorithms[END_REF].

We evaluate the reconstruction results both visually and using the peak signal-to-noise ratio (PSNR) criterion defined for 12 bits dynamic images as

P SN R(x, x) = 20 log 10 4095 1 N x -x 2 , ( 7
)
where N is the number of pixels. We show visual results only for the case of 2.5 bits/pixel as at high coding rate (4.0 bits/pixel) the reconstructed images of the on-board and on-ground chains appear perceptually identical and have an almost equal PSNR. This can be explained by the fact that at this rate, the compression is almost lossless and the coding step could then be "omitted", making both chains equivalent.

The comparison of these two chains in a rate-distortion sense is given figure 4. Visual results are given figures 5 to 8. We can see that for every restoration techniques, an on-board chain always performs better January 9, 2013 DRAFT in term of PSNR. At low coding rate, the difference between the two chains reaches almost 1 dB. Visually, we can check on the figure 5 for example that the on-board chain gives edges which are slightly more blurred than the on-ground chain (particularly visible around the edges of buildings). This is due to the fact that the edges of the image have been enhanced by the deconvolution. The high-frequencies subbands require then more bits to be properly encoded. On smooth area, the on-board chain will reconstructs an image which is more detailled (see figure 8) but on an area which contains many edges, the obtained image will be slightly more deteriorated than the one we would have obtained with an on-ground chain (see figures 5 and 7).

Visually, the FORWARD restoration technique does not give satisfying results and tends to oversmooth the image. If we observe the reconstructed images (figures 5 and 7 for example), we can verify that all the small details are lost. The Stein block thresholding [START_REF] Chesneau | Stein block thresholding for wavelet-based image deconvolution[END_REF] technique seems also to suffer from the same behavior and provides smooth reconstructed images. This method retrieves however more details than the ForWard method but also gives more reconstruction artifacts, see figure 7. The method based on the Shearlets seems to be the most satisfying in term of image quality. It allows to recover the small details of the image without giving too many artifacts.

January 9, 2013 DRAFT From a numerical point of view, the FORWARD restoration algorithm [START_REF] Neelamani | ForWaRD: Fourier-Wavelet Regularized Deconvolution for Ill-Conditioned Systems[END_REF] gives the best PSNR for all rates, though all these methods are slightly outperformed by the restoration technique used by the CNES, at least up to a coding rate of 2.5 bits/pixel. Note that the PSNR of that method is almost constant after the coding rate of 2.5 bits/pixel as the restoration algorithm used by the CNES leaves some residual noise to give the image a physical sense. This residual noise simulates the instrumental noise that one obtains at the outpout of a sensor. This phenomenon only appears from 2.5 bits/pixel, as at this rate the encoder starts to efficiently encode the instrumental noise instead of removing it. Also note that this image characteristic is highly appreciated from image analysis experts. We will use this characteristic to propose a technique to remove the coding artifacts inherent in wavelet-based compression systems. This is the focus of the next part.

(a) (b) (c) (d) (e) (f) (g) (h) (i)

III. QUANTIZING ARTIFACTS REMOVAL

In this part, we briefly describe the state-of-the-art of quantization noise removal methods. Note that this part is only a brief review of quantization noise removal techniques and we will discuss the integration of these techniques in the satellite chain in the last part of the paper.

A. Variational methods for denoising quantization noise

The quantization noise is very specific and visually comes out as structured artifacts which represent the responses of the wavelets. This phenomenon appears as oscillations, basically near the edges of the image (ringing artifacts), when an initial non-null wavelet coefficient has been set to zero by the quantification. Several methods based on a variational framework have been recently proposed [START_REF] Durand | Reconstruction Of Wavelet Coefficients Using Total Variation Minimization[END_REF], [START_REF] Weiss | Compression artifacts reduction using variational methods: algorithms and experimental study[END_REF] to tackle the problem of quantization noise removal. These methods rely on the minimization of an a priori named the Total Variation (TV) [START_REF] Rudin | Nonlinear Total Variation based noise removal algorithms[END_REF]. The TV prior supposes that an image can be modeled as a smooth function with discontinuities along curves. The oscillations created by the artifacts cannot therefore be considered to be natural and do not belong to an image. The simple way to remove these variations is to minimize the l 1 -norm of the gradient of the image, namely the TV, such that these oscillations are replaced by smooth homogeneous regions.

January 9, 2013 DRAFT (a) (b) (c) (d) (e) (f) (g) (h) (i)
Both methods [START_REF] Durand | Reconstruction Of Wavelet Coefficients Using Total Variation Minimization[END_REF] and [START_REF] Weiss | Compression artifacts reduction using variational methods: algorithms and experimental study[END_REF] could globally be formalized as the following minimization problem

x = arg min ∇x 1 subject to x ∈ K , ( 8
)
where x is the denoised image and K is a set that constrains the reconstructed image. The authors of [START_REF] Weiss | Compression artifacts reduction using variational methods: algorithms and experimental study[END_REF] proposed to define the set K such that it constrains the error between the observed and the reconstructed wavelet coefficients to be bounded by the boundaries of the true quantization error.

Let x 0 be the image to encode, W the wavelet transform of the coder and x the decoded image. We denote boundaries of each quantization interval such that

Q
(W x) i = q k , if b k ≤ (W x 0 ) i < b k+1 , ∀i ∈ {0, . . . , N -1}. (9)
We have

b k -q k ≤ (W x 0 ) i -(W x) i < b k+1 -q k , ∀i ∈ {0, . . . , N -1}. (10)
We set α i = b kq k and β i = b k+1q k ; the authors of [START_REF] Weiss | Compression artifacts reduction using variational methods: algorithms and experimental study[END_REF] proposed to define K as the following hypercube

K = x ∈ R N , α i ≤ (W x) i -(W x) i < β i , ∀i ∈ {0, . . . , N -1} . (11)
The true quantization error is unknown as the original image x 0 is not accessible. The bounds α i and β i of this error can however be estimated from the reconstructed image and the knowledge of the quantizing model. So the originality of the method proposed in [START_REF] Weiss | Compression artifacts reduction using variational methods: algorithms and experimental study[END_REF] consists in minimizing the TV of the reconstructed image such that the quantization error belongs to the invervals defined by the boundaries January 9, 2013 DRAFT [START_REF] Kalifa | Deconvolution by thresholding in mirror wavelet bases[END_REF] of the true quantization error. The method proposed in [START_REF] Durand | Reconstruction Of Wavelet Coefficients Using Total Variation Minimization[END_REF] is slightly different and constrains the wavelet coefficients that have not been set to zero by the quantification to remain identical. They define the set K as

K = x ∈ R N , (W x) i = (W x) i , ∀i ∈ M , ( 12
)
where M is the set of coefficients coordinates that have not been set to zero by the quantizing

M = {i ∈ {0, . . . , N -1}, |(W x) i | > 0}. ( 13
)
The idea of the method proposed in [START_REF] Durand | Reconstruction Of Wavelet Coefficients Using Total Variation Minimization[END_REF] is to reconstruct the small coefficients that have been set to zero by the quantification. The method relies on the fact that the minimization of the TV creates flat regions which need small wavelet coefficients to be represented. The presence of the constraint ( 12) is to ensure that only the small coefficients are updated and that the large quantized coefficients, which are likely to be close to the original ones, remain unchanged.

A comparison of the two presented methods is given at the end of this part. We will see however that the flat homogeneous regions created by these methods are not natural in the sense that they cannot be interpreted physically. The problem of quantization noise removal is actually very difficult to address.

The main difficulty lies in the fact that the quantization noise is highly correlated to the signal source and cannot be modeled using classical probability distributions. We present in the next part methods to improve the statistical properties of the quantization noise.

B. Dithering methods for removing quantization artifacts

A dithering technique consists in inserting a noise prior to quantizing to improve the statistics of the quantization error. Among the dithering techniques, we focus on the substractive dithering technique proposed in [START_REF] Schuchman | Dither Signals and Their Effect on Quantization Noise[END_REF] whose particularity is to substract the added noise after quantizing. Let w be the wavelet coefficients of the image to quantize and Q a uniform scalar quantizer of quantizing step ∆. The quantized coefficients w are obtained by applying the quantizing operator Q on the wavelet coefficients w noised by the dithering noise v

w = Q(w + v). ( 14
)
As mentioned previously, a substractive dithering scheme also substracts the dithering noise after the quantizing. The final wavelet coefficients z are

z = w -v = Q(w + v) -v. ( 15
)
January 9, 2013 DRAFT

The authors of [START_REF] Schuchman | Dither Signals and Their Effect on Quantization Noise[END_REF] showed that the global error ǫ = z-w of a substractive dithering system is independent of the system source and is distributed uniformly if the dither noise v can be expressed as a summation of rectangular probability density functions [START_REF] Wannamaker | A theory of nonsubtractive dither[END_REF]. More precisely, the probability density function p ǫ of the global error ǫ can be expressed as [START_REF] Schuchman | Dither Signals and Their Effect on Quantization Noise[END_REF] p

ǫ (ǫ) = ∆Π ∆ (ǫ). [W ∆ * p w * p v ] (-ǫ), (16)
where p v is the probability density function of the dithering noise v, p w is the probability density function of the wavelet coefficients w to quantize and

Π ∆ (ǫ) =        1 ∆ if -∆ 2 < ǫ ≤ ∆ 2 , 0 otherwise (17)
W ∆ (ǫ) = +∞ k=-∞ δ(ǫ -k∆). ( 18
)
The characteristic function (defined as the Fourier transform of the probability density function) is usually easier for studying the correlation properties. It writes [START_REF] Lipshitz | Quantization and Dither: A Theoretical Survey[END_REF]

P ǫ (u) = sinc(u) * [W 1 ∆ (-u)P w (-u)P v (-u)] = +∞ k=-∞ sinc u - k ∆ P w - k ∆ P v - k ∆ = sinc(u) + +∞ k=-∞, k =0 sinc u - k ∆ P w - k ∆ P v - k ∆ . ( 19
)
where

sinc(u) =        sin(π∆u) π∆u , if u = 0 1, otherwise . (20)
From ( 19), we see that the global error ǫ is not function of the source w if the characteristic function P ǫ is reduced to sinc(u). This is verified if the characteristic function of the dither noise P v cancels for

each k ∆ , ∀k ∈ Z * [23] P v k ∆ = 0, ∀k ∈ Z * . ( 21
)
As mentioned by [START_REF] Wannamaker | A theory of nonsubtractive dither[END_REF], the equation ( 21) is satisfied if the dithering noise can be expressed as a summation of rectangular probability density functions [START_REF] Wannamaker | A theory of nonsubtractive dither[END_REF]. In that case, the global error ǫ reduces to an independent and uniformly distributed noise.

This is an encouraging result as it implies that an on-board restoration coupled with a substractive dithering January 9, 2013 DRAFT scheme will result in a restored image with a residual noise which is independent of the original image.

This residual noise can then be interpreted physically as the intrumental noise of the sensor. This aspect of residual noise is very important as it is one of the features seeked from the CNES for the design of restoration methods [START_REF] Dherete | Image de-blurring and application to SPOT5 THR satellite imaging[END_REF]. We will discuss this aspect later as this is the basis of the proposed imaging chain described in section IV.

Finally, we would like to mention also the dithering technique proposed by [START_REF] Stamm | Anti-Forensics of Digital Image Compression[END_REF]. The method proposed in [START_REF] Stamm | Anti-Forensics of Digital Image Compression[END_REF] is slightly different from the classical dithering techniques as it is more focus on the reconstruction of the original wavelet subbands rather than improving the statistics of the quantizing noise. More precisely, the main result of [START_REF] Stamm | Anti-Forensics of Digital Image Compression[END_REF] is that the probability density function of a wavelet subband can be recovered exactly (assuming we know the parameters of its model) from its quantized version by adding a dithering noise v to the quantized coefficients. As previously, we denote by w an original (i.e. prior to quantizing) wavelet subband and w the corresponding quantized subband. The authors of [START_REF] Stamm | Anti-Forensics of Digital Image Compression[END_REF] proposed to model a wavelet subband w (each subband can be treated separately) by a Laplace distribution [START_REF] Li | Text and picture segmentation by the distribution analysis of wavelet coefficients[END_REF]

p w (w) = λ 2 e -λ|w| , (22)
where λ is the scale parameter that can be estimated using classical estimation techniques such that least-squares minimization methods or maximum-likelihood estimations. Similarly to (9), we have

w = q k , if b k ≤ w < b k+1 . (23)
Using the wavelet subband model [START_REF] Rudin | Nonlinear Total Variation based noise removal algorithms[END_REF], we can write that

p w( w = q k ) =                1 2 e -λbk -e -λbk+1 , if k ≥ 1 1 -1 2 e -λb0 -e -λb1 , if k = 0 1 2 e λbk+1 -e λbk , if k ≤ -1. (24)
The output wavelet subband z is given by

z = w + v, ( 25
)
where w is the quantized wavelet subband and v the dithering noise. The output wavelet subband probability density function can be expressed using the law of total probability [START_REF] Stamm | Anti-Forensics of Digital Image Compression[END_REF]

p z (z) = k p z| w(z| w = q k )p w( w = q k ), (26)
January 9, 2013 DRAFT where

p z| w(z| w = q k ) = p v| w(v = z -q k | w = q k ), (27)
is the probability density function of the dither noise v knowing the quantized values w. The authors of [START_REF] Stamm | Anti-Forensics of Digital Image Compression[END_REF] showed that the choices

p v| w(v| w = q k , k = 0) =        1 αk e -sign(qk) λv , if (b k -q k ) ≤ v < (b k+1 -q k ) 0, otherwise (28)
p v| w(v| w = 0) =        1 α0 e -λ|v| , if b 0 > v > b 1 0, otherwise , (29)
with a k being normalization constants and λ an estimated value of the scale parameter λ, lead to the original wavelet subband probability density function, under the condition that λ = λ [27]

p z (z) = k p z| w(z| w = q k )p y ( w = q k ) = k≤-1 1 α k e λ(z-qk) 1 2 e λbk+1 -e λbk 1(b k ≤ z < q k+1 ) + 1 α 0 e -λ|z| 1 - 1 2 e λb0 -e -λb1 1(b 0 ≤ z < b 1 ) + k≥1 1 α k e -λ(z-qk) 1 2 e -λbk -e -λbk+1 1(b k ≤ z < q k+1 ) = λ 2 e -λ|z| = p w (z) (30)
Even if the reconstructed and original subbands will numerically differ, this technique will remove the undesirable observed artifacts, due to the quantization, by filling in the blanks. The fact that we also add dither noise on the null coefficients may also provide the residual noise appreciated from the images expert analysis.

C. Comparison of removal methods for quantization artifacts

We simulate the behavior of the presented quantization removal methods directly on a coded version of the reference (i.e. without any blur or noise) satellite image shown figure 3. The simulation of the complete imaging chain including these techniques is done in the next part. To perform a fair comparison, the image will be coded using the biorthogonal 9/7 wavelet transform [START_REF] Cohen | Biorthogonal bases of compactly supported wavelets[END_REF] followed by the quantizer described in [START_REF] Lipshitz | Quantization and Dither: A Theoretical Survey[END_REF]. As a consequence, the method [START_REF] Stamm | Anti-Forensics of Digital Image Compression[END_REF] has been adapted to this choice. For the substractive dithering method [START_REF] Wannamaker | A theory of nonsubtractive dither[END_REF], we simulated a uniform dithering noise to limit the power of the residual noise.

We only provide visual results as common criteria such as PSNR do not take into account the appreciated physical perception of residual noise.

The results are given figures 9 to 12. Visually, we immediately see that the techniques based on the minimization of the TV create large smooth homogenous regions and remove the small details of the image. This effect is known as the cartoon effect. These flat regions are not considered to be natural for a remote sensing image and are really not appreciated from image analysis experts who clearly prefer a deterioration that can be interpreted physically. As explained previously, this is for example the case of an unstructured residual noise. The substractive dithering system and the method proposed in [START_REF] Stamm | Anti-Forensics of Digital Image Compression[END_REF] give good visual results. Both images are well reconstructed and do not present common artifacts such as ringing or blurry edges. The quality of the image reconstructed with the substractive dithering system actually seems slightly better, particularly on the small details of the image (cars and zebras). As expected, these methods leave a residual noise on the reconstructed image which simulates the instrumental noise that image expert analysis are used to observe.

IV. PROPOSED IMAGING CHAIN

In the previous section, we showed that the dithering techniques may be very interesting to remove the structured artifacts of the coding step. As we have also mentioned in section III, these techniques leave a uniform residual noise which is highly appreciated from the image analysis experts as it can be interpreted physically. More precisely, the "perfect" restored image (as defined by image analysis experts) should owns a residual blur characterized by a target PSF [START_REF] Lambert-Nebout | On-board optical image compression for future high-resolution remote sensing systems[END_REF] along with a uniform residual noise with a fixed standard deviation [START_REF] Dherete | Image de-blurring and application to SPOT5 THR satellite imaging[END_REF].

We presented in section II an on-board restoration technique which gives an image with a residual noise (whose power is very small in comparison to the power of the residual noise obtained from the dithering technique) and a residual blur fully characterized by the target PSF. If we combine these two techniques, i.e. if we use an on-board restoration coupled with a substractive dithering technique, the image obtained at the output of the chain will then present a residual noise (coming from the dithering technique) with the blur of the target PSF (coming from the on-board restoration). And as mentioned previously, a final image with such characteristics is the objective of image analysis experts as it can be interprated as the direct output of an ideal instrument characterized. [START_REF] Durand | Reconstruction Of Wavelet Coefficients Using Total Variation Minimization[END_REF], and (d) is the image obtained using the post-processing technique proposed in [START_REF] Weiss | Compression artifacts reduction using variational methods: algorithms and experimental study[END_REF], (e) is the image reconstructed using the post-processing dithering technique proposed in [START_REF] Stamm | Anti-Forensics of Digital Image Compression[END_REF], (f) is the image reconstructed using the substractive dithering technique [START_REF] Lipshitz | Quantization and Dither: A Theoretical Survey[END_REF] with an uniform dithering noise. The target rate is 2.5 bits/pixel. The image range has been extended to point up the image reconstruction artifacts.

From this remark, we propose the global imaging chain shown figure 13. This chain includes the on-board restoration [START_REF] Patel | Shearlet-Based Deconvolution[END_REF] based on the Shearlets transform and the substractive dithering technique [START_REF] Wannamaker | A theory of nonsubtractive dither[END_REF] to decorrelate the quantizing noise. Note that, in this chain, the quantizer follows the model described in [START_REF] Lipshitz | Quantization and Dither: A Theoretical Survey[END_REF] to respect the substractive dithering scheme hypothesis. The coding step is then decomposed in a 3-levels CDF 9/7 wavelet transform followed by an explicit quantization of the wavelet coefficients and an entropy encoding of the quantized coefficients. The results of this new chain are given figures 14 to 17.

We immediately see that the reconstructed images with the proposed chain do not present any common wavelet compression artifacts (see figures obtained residual noise is however better appreciated from image analysis experts at it can be interpreted physically. More precisely, the proprosed imaging chain produces a reconstructed image which owns the two characteristics of an ideal image: Blur with the target PSF (obtained by the on-board restoration) and a residual unstructured noise [START_REF] Dherete | Image de-blurring and application to SPOT5 THR satellite imaging[END_REF]. The drawback of the proposed method is that the standard deviation of the residual noise is function of the quantizing step while it should be constant for all coding rates.

Consequently, for a low coding rate, the proposed chain gives an image which is more noisy than the instrumental one. Futher works need thus to be done on this aspect.

V. CONCLUSIONS

In this paper, we considered the satellite imaging chain optimization problem. We presented numerical results which showed that the quality of the reconstructed image can be improved if one concedes several changes on the usual design of imaging chain. The first one would be to move the restoration step onboard of the satellite, prior to compression. This allows to reconstruct an image with less reconstruction artifacts, specially on shadows zones. The second deals with the problem of coding noise removal. From the results we presented, we concluded that the current state-of-art algorithms do not give competitve results and that the best option may be to used dithering techniques to transform this stuctured coding noise in an unstructured uniform residual noise. This residual noise is very interesting as it simulates the noise obtained at the output of the instrument, and this property is highly appreciated from photo interpreters. From these conclusions, we propose a new imaging chain and we showed results on a real satellite data. We compared the results of the proposed chain with the ones obtained with the current satellite imaging chain. We showed that the proposed chain gives interesing results and may be particularly efficient at medium and high coding rates (around 3.0 bits/pixel and more). Future works will be focussed to extend our approach such that we get competitive results at low coding rates.

Fig. 1 .

 1 Fig. 1. Processing chain for satellite imaging.

  . Let x be a discrete version of the true analog image, y the instrumental one and x the restored image at the output of the chain. The coding/decoding and restoration are respectively denoted by the operators C and T . We denote by D(x, x) some measure of the distance between the reference image and the restored one and by R(C(y)) some measure of the coding rate of the coded image. The restored image can be expressed as a function of the coding and the restoration by x = T (C(y)). The problem of global optimization consists in finding the optimal C * and T * such that C * , T * = arg min E [D(x, T (C(y)))] subject to C, T R(C(y)) ≤ R c ,

Fig. 2 .

 2 Fig. 2. On-board restoration based satellite imaging chain.

  E[x|y] is the conditional expectation of the original image x knowing the noisy one y. The image E[x|y] can be seen as the best estimator (which minimizes the MSE) of the original image x from y. This image does not depend on the restoration or the compression technique used. So the minimal distortion January 9, 2013 DRAFTD * then writes[START_REF] Wolf | Transmission of noisy information to a noisy receiver with minimum distortion[END_REF]

Fig. 3 .

 3 Fig. 3. Reference image, Cannes harbour (1024 × 1024 pixels).

Fig. 4 .

 4 Fig. 4. Rate-distortion comparison of on-board and on-ground chains in reference to the method currently used by the CNES.

Fig. 5 .Fig. 6 .

 56 Fig. 5. Visual comparison of on-board and on-ground chains. (a) is the reference image, (b) is the instrumental image (output of the acquisition, P SN R = 32.69 dB), (c) is the reconstructed image provided by the CNES (P SN R = 45.93 dB), (d) and (e) are the reconstructed images respectively from the Shearlets based on-board (P SN R = 45.67 dB) and on-ground (P SN R = 45.12 dB) chains, (f) and (g) are the reconstructed images respectively from the block thresholding based on-board (P SN R = 45.38 dB) and on-ground (P SN R = 44.76 dB) chains, (h) and (i) are the reconstructed images respectively from the FORWARD based on-board (P SN R = 45.83 dB) and on-ground (P SN R = 45.28 dB) chains. The target rate is 2.5 bits/pixel. The image range has been extended to point up the image reconstruction artifacts.

Fig. 7 .

 7 Fig. 7. Visual comparison of on-board and on-ground chains. (a) is the reference image, (b) is the instrumental image (output of the acquisition, P SN R = 32.69 dB), (c) is the reconstructed image provided by the CNES (P SN R = 45.93 dB), (d) and (e) are the reconstructed images respectively from the Shearlets based on-board (P SN R = 45.67 dB) and on-ground (P SN R = 45.12 dB) chains, (f) and (g) are the reconstructed images respectively from the block thresholding based on-board (P SN R = 45.38 dB) and on-ground (P SN R = 44.76 dB) chains, (h) and (i) are the reconstructed images respectively from the FORWARD based on-board (P SN R = 45.83 dB) and on-ground (P SN R = 45.28 dB) chains. The target rate is 2.5 bits/pixel. The image range has been extended to point up the image reconstruction artifacts.

Fig. 8 .

 8 Fig. 8. Visual comparison of on-board and on-ground chains. (a) is the reference image, (b) is the instrumental image (output of the acquisition, P SN R = 32.69 dB), (c) is the reconstructed image provided by the CNES (P SN R = 45.93 dB), (d) and (e) are the reconstructed images respectively from the Shearlets based on-board (P SN R = 45.67 dB) and on-ground (P SN R = 45.12 dB) chains, (f) and (g) are the reconstructed images respectively from the block thresholding based on-board (P SN R = 45.38 dB) and on-ground (P SN R = 44.76 dB) chains, (h) and (i) are the reconstructed images respectively from the FORWARD based on-board (P SN R = 45.83 dB) and on-ground (P SN R = 45.28 dB) chains. The target rate is 2.5 bits/pixel. The image range has been extended to point up the image reconstruction artifacts.

Fig. 9 .

 9 Fig. 9. Visual comparison of quantizing removal techniques. (a) is the reference image, (b) is the compressed image, (c) is the image obtained using the post-processing technique proposed in [20], and (d) is the image obtained using the post-processing

Fig. 10 .

 10 Fig. 10. Visual comparison of quantizing removal techniques. (a) is the reference image, (b) is the compressed image, (c) is the image obtained using the post-processing technique proposed in[START_REF] Durand | Reconstruction Of Wavelet Coefficients Using Total Variation Minimization[END_REF], and (d) is the image obtained using the post-processing technique proposed in[START_REF] Weiss | Compression artifacts reduction using variational methods: algorithms and experimental study[END_REF], (e) is the image reconstructed using the post-processing dithering technique proposed in[START_REF] Stamm | Anti-Forensics of Digital Image Compression[END_REF], (f) is the image reconstructed using the substractive dithering technique[START_REF] Lipshitz | Quantization and Dither: A Theoretical Survey[END_REF] with an uniform dithering noise. The target rate is 2.5 bits/pixel. The image range has been extended to point up the image reconstruction artifacts.

Fig. 11 .

 11 Fig. 11. Visual comparison of quantizing removal techniques. (a) is the reference image, (b) is the compressed image, (c) is the image obtained using the post-processing technique proposed in[START_REF] Durand | Reconstruction Of Wavelet Coefficients Using Total Variation Minimization[END_REF], and (d) is the image obtained using the post-processing technique proposed in[START_REF] Weiss | Compression artifacts reduction using variational methods: algorithms and experimental study[END_REF], (e) is the image reconstructed using the post-processing dithering technique proposed in[START_REF] Stamm | Anti-Forensics of Digital Image Compression[END_REF], (f) is the image reconstructed using the substractive dithering technique[START_REF] Lipshitz | Quantization and Dither: A Theoretical Survey[END_REF] with an uniform dithering noise. The target rate is 2.5 bits/pixel. The image range has been extended to point up the image reconstruction artifacts.

Fig. 12 .

 12 Fig. 12. Visual comparison of quantizing removal techniques. (a) is the reference image, (b) is the compressed image, (c) is the image obtained using the post-processing technique proposed in[START_REF] Durand | Reconstruction Of Wavelet Coefficients Using Total Variation Minimization[END_REF], and (d) is the image obtained using the post-processing technique proposed in[START_REF] Weiss | Compression artifacts reduction using variational methods: algorithms and experimental study[END_REF], (e) is the image reconstructed using the post-processing dithering technique proposed in[START_REF] Stamm | Anti-Forensics of Digital Image Compression[END_REF], (f) is the image reconstructed using the substractive dithering technique[START_REF] Lipshitz | Quantization and Dither: A Theoretical Survey[END_REF] with an uniform dithering noise. The target rate is 2.5 bits/pixel. The image range has been extended to point up the image reconstruction artifacts.

Fig. 13 .Fig. 14 .Fig. 15 .Fig. 16 .Fig. 17 .

 1314151617 Fig. 13. Proposed satellite imaging chain

  The second problem is to design a coder C which minimizes the coding error.

	of global joint optimization (2) is to use an on-board restoration such that the global optimization problem
	can be splitted in two independent ones. The first problem is to optimize the on-board restoration such
	that it is close to E[x|y].
	January 9, 2013	DRAFT

January 9, 2013 DRAFT