Caroline Chaux
email: chaux@univ-mlv.fr

Jean-Christophe Pesquet
email: pesquet@univ-mlv.fr

Laurent Duval
email: laurent.duval@ifp.fr

2D dual-tree complex biorthogonal M-band wavelet transform

Keywords: Wavelet transforms, Hilbert transforms, Image analysis, Image processing, Gaussian noise

published or not. The documents may come   L'archive ouverte pluridisciplinaire

INTRODUCTION

Natural images often require a transformed domain processing in order to improve their quality (denoising or restoration tasks) or to be efficiently compressed. Powerful tools allowing improved representations of images have been developed in the last ten years. They aim at capturing directional features and often introduce some redundant analysis. They can be classified into different categories: steerable filters, "lets" transforms including bandelets [START_REF] Pennec | Sparse geometric image representations with bandelets[END_REF], curvelets [START_REF] Candès | Curves and Surfaces, chapter Curvelets -a surprisingly effective nonadaptive representation for objects with edges[END_REF], contourlets [START_REF] Do | The contourlet transform: an efficient directional multiresolution image representation[END_REF], ... Meanwhile, the dual-tree wavelet decomposition has been introduced by N. Kingsbury [5]. This transform is based on a combination of classical wavelet decompositions. It has been further investigated by I. Selesnick [START_REF] Selesnick | Hilbert transform pairs of wavelet bases[END_REF]. The standard real dual-tree decomposition is 2 times redundant and is nearly shift-invariant. In the 2-band case, in two dimensions, a complex version of this decomposition has been obtained (then, the redundancy becomes equal to 4) and, more recently, it has also been extended to biorthogonal representations [START_REF] Yu | Hilbert transform pairs of biorthogonal wavelet bases[END_REF]. The dual-tree transform is interesting for several reasons: good directional analysis, low redundancy, improved shift-invariance property, simplicity of implementation, reduced computational cost,... However, a limiting factor lies in its dyadic structure which introduces some inherent design constraints. For example, it is well-known that is impossible to obtain real, orthogonal, symmetric and compactly supported dyadic wavelets, except for the Haar one. In one of our recent works [START_REF] Chaux | Image analysis using a dual-tree M -band wavelet transform[END_REF], we have extended the real dual-tree transform to the M -band case, which allowed us to gain more freedom in the choice of the filters while providing a more accurate frequency analysis.

In this work, we propose to extend the real 2D dual-tree M -band wavelet transform to the biorthogonal complex case, for an arbitrary integer value of the decimation factor M . The paper is organized as follow: in Section 2 we show how to design Hilbert pairs of biorthogonal wavelets. Then, Section 3 is devoted to the complex decomposition. Finally, denoising applications are reported in Section 4.

Throughout the paper, the following notations will be used: let M be an integer greater than or equal to 2, N M = {0, . . . , M -1} and N ⋆ M = {1, . . . , M -1}. Besides, a denotes the Fourier transform of a function a, ⌈u⌉ denotes the upper integer part of a real u, (δ m ) m∈Z is the Kronecker sequence (equal to 1 if m = 0 and 0 otherwise).

BIORTHOGONAL HILBERT PAIRS OF WAVELETS

We address in this section the problem of constructing a Hilbert pair of M -band biorthogonal wavelets in the one dimensional case. A biorthogonal wavelet basis corresponding to a multiresolution analysis of L 2 (R) is associated to an analysis and a synthesis filter bank and the problem is then equivalent to the design of 2M filters.

Problem statement

Let us start from an M -band biorthogonal wavelet decomposition of L 2 (R). This decomposition is based on the joint use of two sets of basis functions: (ψ m ) 0≤m<M and ( ψ m ) 0≤m<M which satisfy the following scaling equations expressed in the frequency domain:

∀m ∈ N M , √ M ψ m (M ω) = H m (ω) ψ 0 (ω), √ M ψ m (M ω) = H m (ω) ψ 0 (ω). (1)
Furthermore, the associated M -band filter banks with frequency responses (H m ) 0≤m<M and ( H m ) 0≤m<M satisfy per-fect reconstruction properties:

∀(m, m ′ ) ∈ N 2 M , M -1 p=0 H m (ω + p 2π M )H * m ′ (ω + p 2π M ) = M δ m-m ′ . (2)
Similarly to the orthogonal case [START_REF] Chaux | Image analysis using a dual-tree M -band wavelet transform[END_REF], we want to construct Hilbert pairs of wavelets by defining functions

(ψ H m ) 0≤m<M and ( ψ H m ) 0≤m<M such that ∀m ∈ N ⋆ M , ψ H m (ω) = -ı sign(ω) ψ m (ω) ψ H m (ω) = -ı sign(ω) ψ m (ω), (3)
where sign designates the signum function. The filters with frequency responses (G m ) 0≤m<M and ( G m ) 0≤m<M associated with these dual basis functions must also satisfy perfect reconstruction conditions:

∀(m, m ′ ) ∈ N 2 M , M -1 p=0 G m (ω + p 2π M )G * m ′ (ω + p 2π M ) = M δ m-m ′ . ( 4
)
We now study how to design such filters.

Dual filter construction

The Hilbert conditions (3) imply that:

∀m ∈ N ⋆ M , | ψ H m (ω)| = | ψ m (ω)| and | ψ H m (ω)| = | ψ m (ω)|.
In

∀m ∈ N M , G m (ω) = e -ıθm(ω) H m (ω) G m (ω) = e -ı θm(ω) H m (ω), (5)
where θ m and θ m are real-valued 2π-periodic phase functions. Assuming that ( 2) is satisfied, one can easily check that the perfect reconstruction conditions (4) for the dual filter bank are fulfilled if for all (m, m

′ ) ∈ N 2 M , θ m = θ m (mod 2π) and θ m ′ -θ m = θ m ′ -θ m is 2π/M -periodic.
Conditions to be satisfied by θ m and θ m are thus similar to those obtained in the orthogonal case [START_REF] Chaux | Image analysis using a dual-tree M -band wavelet transform[END_REF]. It can be deduced that

∀p ∈ 0, . . . , M 2 -1 , ∀ω ∈ p 2π M , (p + 1) 2π M , θ 0 (ω) = θ 0 (ω) = (d + 1 2 )(M -1) ω -pπ, and ∀m ∈ {1, . . . , M -1}, θ m (ω) = θ m (ω) =    π 2 -d + 1 2 ω if ω ∈]0, 2π[, 0 if ω = 0,
where d ∈ Z.

COMPLEX 2D DUAL-TREE WAVELET TRANSFORM

Our objective in this section is to extend the complex dual tree transform to the M -band case. We will see that an increased number of directions can be selected and show how to implement the associated decomposition. Finally, the problem of finding an optimized reconstruction will be discussed.

Direction selection

For all m ∈ N ⋆ M and ǫ ∈ {-1, 1}, let us define the following complex-valued functions from the wavelets (ψ m ) m∈N ⋆ M used for the analysis stage:

∀x ∈ R, ψ ǫ m (x) = ψ m (x) + ıǫψ H m (x).
If ǫ = 1 (resp. ǫ = -1), this corresponds to an analytic (resp. anti-analytic) wavelet. By considering all the possible combinations of tensor products of such wavelets, we obtain:

for all (m 1 , m 2 ) ∈ N ⋆2 M , (ǫ 1 , ǫ 2 ) ∈ {-1, 1} 2 and (x 1 , x 2 ) ∈ R 2 , ψ ǫ1 m1 (x 1 )ψ ǫ2 m2 (x 2 ) = ψ m1 (x 1 )ψ m2 (x 2 ) -ǫ 1 ǫ 2 ψ H m1 (x 1 )ψ H m2 (x 2 ) + ı ǫ 1 ψ H m1 (x 1 )ψ m2 (x 2 ) + ǫ 2 ψ m1 (x 1 )ψ H m2 (x 2 ) (6)
which reads in the Fourier domain:

ψ ǫ1 m1 (ω 1 ) ψ ǫ2 m2 (ω 2 ) = (1 + ǫ 1 sign(ω 1 ))(1 + ǫ 2 sign(ω 2 )) ψ m1 (ω 1 ) ψ m2 (ω 2 ).
The latter expression shows that, depending on the choice of ǫ 1 and ǫ 2 , it is possible to select each quadrant of the frequency plane as represented in Fig. 1.

Let us consider the real part (RP) and the imaginary part (IP) of the expression in [START_REF] Selesnick | Hilbert transform pairs of wavelet bases[END_REF] and consider once again all possible situations. This leads to the two following cases:

➢ ǫ 1 = ǫ 2 = 1 or ǫ 1 = ǫ 2 = -1 : RP: ψ m1 (x 1 )ψ m2 (x 2 ) -ψ H m1 (x 1 )ψ H m2 (x 2 ) IP: sign(ǫ 1 ) ψ H m1 (x 1 )ψ m2 (x 2 ) + ψ H m2 (x 2 )ψ m1 (x 1 ) 00
➢ ǫ 1 = -ǫ 2 = 1 or ǫ 1 = -ǫ 2 = -1: RP: ψ m1 (x 1 )ψ m2 (x 2 ) + ψ H m1 (x 2 )ψ H m2 (x 2 ) IP: sign(ǫ 1 )(ψ H m1 (x 1 )ψ m2 (x 2 ) -ψ H m2 (x 2 )ψ m1 (x 1 )).
In other words, in order to implement the decomposition onto (ψ ǫ1 m1 (x 1 )ψ ǫ2 m2 (x 2 )) (m1,m2)∈N ⋆2 M , we need four separable 2D wavelet analyses (see the corresponding filter bank structures in Fig. 2). The two upper analyses, which are also present in the real parts, are implemented in the real biorthogonal case.

These four M -band wavelet decompositions generate coefficients denoted by c

•• j,m [k], c HH j,m [k], c H• j,m [k] and c •H j,m [k] where j ∈ Z denotes the resolution level, m = (m 1 , m 2 ) ∈ N 2
M is the frequency band index and k ∈ Z 2 is the spatial position. Finally, to obtain the desired directional analysis, we compute the inner product of the image to be processed with the RP and IP terms of the two-dimensional wavelets as previously expressed. This is equivalent to perform the following linear combinations of the subbands: for all m ∈ N ⋆2 M ,

d •• j,m [k] = 1 √ 2 (c •• j,m [k] + c HH j,m [k]) d HH j,m [k] = 1 √ 2 (c •• j,m [k] -c HH j,m [k]) d H• j,m [k] = 1 √ 2 (c H• j,m [k] + c •H j,m [k]) d •H j,m [k] = 1 √ 2 (c H• j,m [k] -c •H j,m [k]).

Prefiltering stage

In a digital implementation of the dual-tree decomposition, prefilters must be added to make the transition from the analog formalism to discrete processing [START_REF] Chaux | Image analysis using a dual-tree M -band wavelet transform[END_REF].

H * M -1 H * M -1 G * 0 H * 0 ↓ (M, M ) H 0 H 1 H 0 H 0 H 0 G 0 H 0 G 1 G * 0 H * 1 G * 0 G * 1 H * M -1 G * M -1 G * M -1 H * M -1 G * M -1 G * M -1 ↓ (M, M ) ↓ (M, M ) ↓ (M, M ) ↓ (M, M ) ↓ (M, M ) ↓ (M, M ) ↓ (M, M ) ↓ (M, M ) ↓ (M, M ) ↓ (M, M ) ↓ (M, M ) ↑ (M, M ) ↑ (M, M ) ↑ (M, M ) ↑ (M, M ) ↑ (M, M ) ↑ (M, M ) ↑ (M, M ) ↑ (M, M ) ↑ (M, M ) ↑ (M, M ) ↑ (M, M ) ↑ (M, M ) H M -1 H M -1 G M -1 G M -1 G M -1 H M -1 H M -1 G M -1 H * 0 H * 0 G 0 G 0 G 0 G 1 H * 0 G * 0 H * 0 G * 1 G 0 H 0 G 0 H 1 H * 0 H * 1 G * 0 G * 0 Fig. 2.
Two-dimensional analysis/synthesis filter banks corresponding to a biorthogonal complex wavelet transform.

Assume that the digital image (r k1,k2 ) (k1,k2) to be analyzed is related to its analog version r(x 1 , x 2 ) by the relation:

r(x 1 , x 2 ) = (k1,k2) r k1,k2 χ(x 1 -k 1 , x 2 -k 2 )
where χ is an interpolation function in L 2 (R 2 ). To perform the first separable wavelet decomposition, we need to determine the approximation coefficients: c 0,(0,0

) [k 1 , k 2 ] = r(x 1 , x 2 ), ψ 0 (x 1 -k 1 )ψ 0 (x 2 -k 2 )
. Similar approximations have to be computed for the three remaining trees. It is easy to show that these operations can be performed by applying four digital prefilters with frequency responses:

F 1 (ω 1 , ω 2 ) = ∞ p1=-∞ ∞ p2=-∞ χ(ω 1 + 2p 1 π, ω 2 + 2p 2 π) ψ * 0 (ω 1 + 2p 1 π) ψ * 0 (ω 2 + 2p 2 π) F 2 (ω 1 , ω 2 ) = ∞ p1=-∞ ∞ p2=-∞ χ(ω 1 + 2p 1 π, ω 2 + 2p 2 π) ( ψ H 0 (ω 1 + 2p 1 π)) * ( ψ H 0 (ω 2 + 2p 2 π)) * F 3 (ω 1 , ω 2 ) = ∞ p1=-∞ ∞ p2=-∞ χ(ω 1 + 2p 1 π, ω 2 + 2p 2 π) ( ψ H 0 (ω 1 + 2p 1 π)) * ψ * 0 (ω 2 + 2p 2 π) F 4 (ω 1 , ω 2 ) = ∞ p 1 =-∞ ∞ p 2 =-∞ χ(ω 1 + 2p 1 π, ω 2 + 2p 2 π) ψ * 0 (ω 1 + 2p 1 π)( ψ H 0 (ω 2 + 2p 2 π)) * .

Reconstruction

As already mentioned, the complex dual-tree transform has a redundancy of a factor 4. As a consequence, the reconstruction is not unique and we have to pay attention to the choice of the synthesis scheme. Let r ∈ ℓ 2 (Z 2 ) be the vector of image values and c •• , c HH , c H• , c •H denote the coefficient vectors generated by the complex dual-tree transform. The linear combination of the subbands is not taken into account as it reduces to a basic isometry. The global decomposition operator can be written as

D c : r →     c •• c HH c H• c •H     =     D 1 r D 2 r D 3 r D 4 r     (7)
where, for all i ∈ {1, . . . , 4}, D i = U i F i where F i represents the i-th prefiltering operation and U i is the i-th M -band separable biorthogonal wavelet decomposition.

A robust reconstruction of r is then obtained as the solution of the minimization problem

inf r 4 i=1 c i -D i r 2 Q i (8)
where Q i is a positive self-adjoint operator and . 2 Qi = ., Q i . . The minimizer allows us to defines a generalized pseudo-inverse of the operator D c which is expressed as:

D c ♯ = 4 i=1 F † i U † i Q i U i F i -1 (F † 1 U † 1 Q 1 , . . . , F † 4 U † 4 Q 4 ).
In practice, iterative approaches are generally necessary to compute D c ♯ . Notice however that, by choosing

Q i = (U i U † i ) -1
, the generalized pseudo-inverse takes a simple form, which can be easily implemented by combining the standard synthesis filter bank structures with filtering operations.

EXPERIMENTAL RESULTS

We apply these different decompositions to image denoising. More precisely, we aim at restoring the 512 × 512 Barbara image which is corrupted by an additive zero-mean white Gaussian noise. Four orthogonal decompositions are tested: the real dual-tree transform with dyadic (DTT2) and 4-band (DTT4) filter banks; the complex dual-tree transform with dyadic (DTT2cx) and 4-band (DTT4cx) filter banks. In this case, we use symlets of length 8 and 4-band filters of length 22. Moreover, we consider 4 additional biorthogonal transforms: the dual-tree transform with dyadic (DTTbi2) and 4band (DTTbi4) filter banks; the complex dual-tree transform with dyadic (DTTbi2cx) and 4-band (DTTbi4cx) filter banks. In this case, we use standard biorthogonal 5/3 and 9/7 wavelets and the associated 4-band basis functions corresponding to an equal-subband wavelet packet analysis. Classical estimators are applied to the resulting coefficients: Visushrink (Visu), SUREshrink (SURE) and the bivariate shrinkage (Biv) [START_REF] Selesnick | Bivariate shrinkage with local variance estimation[END_REF]. Quantitative results are given in Table 1. Here, the initial SNR is equal to 5.67dB. One can observe that, for all estimators, the 4-band structures always bring significant improvements (at least 0.2dB) with respect to the dyadic ones. Moreover, the complex decompositions outperform the real ones; the differences are however smaller with the bivariate shrinkage, for 4-band orthogonal filter banks. One can also see that the 9/7 biorthogonal decompositions provide results close to orthogonal ones. Nevertheless, biorthogonal decompositions could be more useful in image coding applications where they are known to be very effective.

Visu

  addition, by imposing that | ψ H 0 (ω)| = | ψ 0 (ω)|, and | ψ H 0 (ω)| = | ψ 0 (ω)|, the scaling equations (1) and their counterparts for the dual wavelets lead to

2 Fig. 1 .

 21 Fig. 1. Direction selection in the 2D frequency plane for M = 4, j = 1 and (m 1 , m 2 ) = (2, 1). The four crosshatched areas are separated using tensor products of analytic/anti-analytic wavelets.

Table 1 .

 1 Denoising results in terms of SNR (in dB).

				SURE	Biv
		DTT2	9.47	12.66 13.69
		DTT4	10.70 13.23 14.31
		DTT2cx	9.72	12.84 13.95
		DTT4cx	11.06 13.41 14.32
		DTTbi2	8.55	11.65 12.73
	5/3	DTTbi4	9.39	12.14 13.45
		DTTbi2cx	8.78	11.88 13.44
		DTTbi4cx	9.71	12.35 13.82
		DTTbi2	9.33	12.40 13.62
	9/7	DTTbi4	10.55 12.96 14.25
		DTTbi2cx	9.56	12.59 13.97
		DTTbi4cx 10.90 13.13 14.34


Caroline Chaux

Patrick Louis Combettes

Jean-Christophe Pesquet

Valérie R Wajs

Iterative image deconvolution using overcomplete representations

come


Caroline Chaux

Laurent Duval

Amel Benazza-Benyahia

Jean-Christophe Pesquet

A new estimator


Caroline Chaux
email: chaux@univ-mlv.fr

Laurent Duval
email: laurent.duval@ifp.fr

Jean-Christophe Pesquet
email: pesquet@univ-mlv.fr

HILBERT PAIRS OF M-BAND ORTHONORMAL WAVELET BASES

. We also show that, when the primal filter bank and its wavelets have symmetry, it is inherited by their duals. Furthermore, we give a design example where the number of vanishing moments of the approximate dual wavelets is imposed numerically to be the same as for the primal ones.

INTRODUCTION

Critically subsampled two-channel filter banks with perfect reconstruction are proved to be efficient tools in the field of signal or higher dimensional data processing, especially for denoising, detection, or compression. Wavelets have demonstrated excellent performance in numerous applications but their successes are mitigated by some natural limitations: a lack of shift-invariance and some design restrictions.

On the one hand, a classical drawback of the real discrete wavelet transform (DWT) is its shift-variance and its poor directionality in higher dimensions. One possible solution comes from the adjunction of some redundancy in the transform. Shift-invariant wavelet transforms [START_REF] Pesquet | Time invariant orthonormal representations[END_REF] suppress shift dependencies, at the expense of an increased computational cost and a redundancy factor depending on the decomposition level, typically log 2 N in one dimension, N being the number of samples. Less demanding solutions have been developed, by the use of complex filters for real signals (we refer to [START_REF] Zhang | Orthogonal complex filter banks and wavelets: some properties and design[END_REF] for an overview and design examples), or by employing other wavelet frames. One of the most interesting and effective proposal is the dual-tree discrete wavelet transform, proposed by N. Kingbury [START_REF] Kingsbury | The dual-tree complex wavelet transform: a new technique for shift invariance and directional filters[END_REF]. In the latter and companion works, two wavelet trees are developed in parallel, with root filters being offset by one half sample to each other. The resulting analysis yields a limited redundancy of 2 d for d-dimensional signals, with a much lower shift sensitivity than the DWT and much better directionality. The design of the associated filters was recently addressed and extended by I. Selesnick [START_REF] Selesnick | Hilbert transform pairs of wavelet bases[END_REF], through an -at least approximate -Hilbert pair formulation for the "dual" wavelets. The later result was recently proved to be the "only" possible solution in [START_REF] Özkaramanli | On the phase condition and its solution for Hilbert transform pairs of wavelets bases[END_REF]. Selesnick also proposed the double-density DWT, with a redundancy factor of (3 d -1)/(2 d -1). He finally combined both frame approaches in [START_REF] Selesnick | The double-density dual-tree discrete wavelet transform[END_REF]. An extension of the dual-tree DWT named phaselets has recently been introduced by R. Gopinath in [START_REF] Gopinath | The phaselet transform -an integral redundancy nearly shift-invariant wavelet transform[END_REF]. The advantages of Hilbert pairs have been earlier recognized by other authors, such as P. Abry et al. [START_REF] Abry | Multiresolution transient detection[END_REF] or W. Freeman et al. [START_REF] Freeman | The design and use of steerable filters[END_REF]. Wavelet frames are also addressed by oversampling, as in the context of shiftable transforms [START_REF] Simoncelli | Shiftable multi-scale transforms[END_REF]. More recently, several authors have proposed a projection scheme with an explicit control of the redundancy [START_REF] Fernandes | Complex wavelet transforms with allpass filters[END_REF]. The developed framework resulted in the original 3-band complex filter bank described in [START_REF] Van Spaendonck | Orthogonal Hilbert transform filter banks and wavelets[END_REF]. Another recent approach is the "analytic" wavelet thresholding developped by S. Ohlede et al. [START_REF] Olhede | Analytic wavelet thresholding[END_REF], in the context of denoising.

On the other hand, classical 2-band wavelets also suffer from constraints in their design. In image applications, filter symmetry is often considered as important to minimize boundary artifacts and save computations, as well as finite support for dealing with realistic finite length signals. But compactly supported, real, orthogonal wavelets cannot have linear phase, except for the Haar system, whose regularity is far too limited. As a consequence, some authors have developed quickly decaying wavelets, nearly symmetric filters (such as the symmlets), symmetric biorthogonal filters (such as the 9/7 or 5/3 used in the JPEG 2000 standard) or "almost real" complex filters (i.e. with relatively few energy remaining in the imaginary component of the filtered signal, e.g. in [START_REF] Belzer | Complex, linear-phase filters for efficient image coding[END_REF]). Although having been very successful in many applications, the dyadic decomposition is sometimes considered relatively coarse in the highest frequency subbands. Furthermore, since the high-pass filter in the 2-band scheme heavily depends on the low-pass filter, the wavelet filters' design sometimes lacks degrees of freedom. These constraints have altogether motivated the development of more general M-band structures such as M-band wavelets [START_REF] Steffen | Theory of regular M-band wavelet bases[END_REF] or M-channel Lapped Transforms.

Based on the aforementioned shortcomings and potential workouts to alleviate some of the DWT drawbacks, we propose in this work to extend the dual-tree framework to Mband wavelets. In other words, we develop a redundant decomposition on two M-band orthogonal wavelet bases, where the M -1 wavelets are required to be Hilbert pairs. The organization of the paper is as follow : in Section 2, we explain how to construct the related M-band filter bank. In Section 3, we provide a design example using the resulting conditions. Some conclusions are drawn in Section 4.

G * 0 G * 1 G * M-1 ↓ M ↓ M ↓ M ↑ M ↑ M ↑ M G 0 G 1 G M-1 H * 0 H * 1 H * M-1 ↓ M ↓ M ↓ M ↑ M ↑ M ↑ M H 0 H 1 H M-1
Figure 1: A pair of analysis/synthesis M-band para-unitary filter banks.

CONSTRUCTION OF M-BAND HILBERT PAIRS

Problem statement

Let M be an integer greater than or equal to 2. Recall that an M-band multiresolution analysis is defined by one scaling function ψ 0 ∈ L 2 (R) and (M -1) mother wavelets ψ m ∈ L 2 (R), m ∈ {1, . . . , M -1} [15]. These functions are solutions of the following scaling equations:

∀m ∈ {0, . . . , M -1}, 1 √ M ψ m t M = ∞ ∑ k=-∞ h m [k]ψ 0 (t -k), (1)
where the sequences (h m [k]) k∈Z are real-valued and square summable. The Fourier transform of (h m [k]) k∈Z is a 2πperiodic function, denoted by H m . For the set of functions ∪ M-1 m=1 {M -j/2 ψ m (M -j t -k), ( j, k) ∈ Z 2 } to correspond to an orthonormal basis of L 2 (R), the following para-unitarity conditions must hold:

∀(m, m ′ ) ∈ {0, . . . , M -1} 2 , M-1 ∑ p=0 H m ω + p 2π M H * m ′ ω + p 2π M = Mδ m-m ′ , (2)
where δ m = 1 if m = 0 and 0 otherwise. In this case, cas- cading the M-band para-unitary analysis and synthesis filter bank, depicted in Fig. 1, allows to decompose and to reconstruct perfectly a given signal.

Our objective is to construct a "dual" M-band multiresolution analysis defined by a scaling function ψ H 0 and mother wavelets ψ H m , m ∈ {1, . . . , M -1}. More precisely, the mother wavelets will be obtained by a Hilbert transform from the "original" wavelets ψ m , m ∈ {1, . . . , M -1}. In the Fourier domain, the desired property reads:

∀m ∈ {1, . . . , M -1}, ψ H m (ω) = -ı sign(ω) ψ m (ω), (3
) where a denotes the Fourier transform of a function a and sign is the signum function. Furthermore, the functions ψ H m are defined by scaling equations similar to (1) involving realvalued sequences (g m [k]) k∈Z . In order to generate a dual M-band orthonormal wavelet basis of L 2 (R), the Fourier transforms G m of these sequences must also satisfy the paraunitarity conditions as expressed by Eq. (2).

Sufficient conditions

If we further impose that:

| ψ H 0 (ω)| = | ψ 0 (ω)|, (4)
the scaling equations combined with the Hilbert condition (3) lead to

∀m ∈ {0, . . . , M -1}, G m (ω) = e -ıθ m (ω) H m (ω), (5)
where θ m is 2π-periodic. These functions should also be odd (for real filters) and thus only need to be determined over [0, π]. Assuming that ( 2) is satisfied, it is then straightforward to verify that the para-unitarity conditions for the dual filter bank hold if the following condition holds:

(P). for all (m, m ′ ) ∈ {0, . . . , M -1} 2 with m < m ′ , α m,m ′ = θ m ′ -θ m is a 2π/M-periodic function.
We are then able to state the following result whose proof is omitted here:

Proposition 1. Assume that Conditions (5) hold. A necessary and sufficient condition for (3) to be satisfied is:

∀m ∈ {1, . . . , M -1}, α 0,m ω M + β (ω) = π 2 sign(ω)[2π], (6)
where

β (ω) = ∞ ∑ i=1 θ 0 ω M i . (7)
The above two relations constitute a generalization to the M-band case of a famous result by Selesnick [START_REF] Selesnick | Hilbert transform pairs of wavelet bases[END_REF] restricted to dyadic wavelets. It is worth noting that the function β is given by the following "additive" scaling equation:

β (ω) = β ω M + θ 0 ω M . ( 8
)

Linear phase solution

In the 2-band case (under weak assumptions), θ 0 verifying Eqs. ( 6) and ( 7) is a linear function on [-π, π[ [START_REF] Selesnick | Hilbert transform pairs of wavelet bases[END_REF]. In the Mband case, we will slightly restrict this constraint on a smaller interval by imposing:

∀ω ∈ [0, 2π/M[, θ 0 (ω) = γω, (9)
where γ ∈ R. Using Eq. ( 8), after some tedious calculations, we prove the following result.

Proposition 2. Under the three conditions (5), (P) and (9), there exists a unique solution to Eq. (6) (mod. 2π), given by

∀m ∈ {1, . . . , M -1}, ∀ω ∈]0, 2π/M[, α 0,m (ω) = π 2 - M 2 ω, ( 10
)
and ∀p ∈ 0, . . . , M 2 -1 , ∀ω ∈ p 2π M , (p + 1) 2π M , θ 0 (ω) = M -1 2 ω -pπ. ( 11
)
It should be noted that unlike in the 2-band case, θ 0 may exhibit discontinuities on ]0, π[ due to the pπ term. We subsequently deduce that the functions θ m , m ∈ {1, . . . , M -1}, should be chosen all equal to

θ m (ω) = π 2 - ω 2 if ω ∈]0, 2π[, 0 if ω = 0. ( 12
)

Symmetry properties

As already pointed out, one of the main advantage of the M-band case with M > 2 is to allow the construction of non-trivial real orthonormal bases with compact support and symmetric (or antisymmetric) wavelets. Assume that these properties are fulfilled for the primal filter bank. Although the Hilbert transformed wavelets are no longer compactly supported, we now show that the dual filters and wavelets still inherit these symmetry properties. Indeed, it can be proved that: We deduce from this result that, if ψ 0 is symmetric w.r.t.

Proposition 3.
τ = L -1 2(M -1) (13)
and, for m ∈ {1, . . . , M -1}, ψ m is symmetric (resp. antisym- metric) w.r.t. τ, then ψ H 0 is symmetric w.r.t. τ + 1/2 and, for m ∈ {1, . . . , M -1}, ψ H m is antisymmetric (resp. symmetric) w.r.t. τ.

DESIGN EXAMPLE

To illustrate the interest of the above conditions, we consider a 4-band (M = 4) para-unitary filter bank and design its dual filter bank.

Design of the Filters

The primal filter bank is taken from [START_REF] Alkin | Design of efficient M-band coders with linearphase and perfect-reconstruction properties[END_REF]. It consists in four finite impulse response (FIR) filters of length 16, generating regular, orthonormal and symmetric basis functions. The corresponding filter coefficients are listed in the first four columns of Table 1.

The frequency responses of these filters can be calculated and the frequency responses of the dual filter bank can be deduced using Eqs. ( 5), [START_REF] Fernandes | Complex wavelet transforms with allpass filters[END_REF] and [START_REF] Van Spaendonck | Orthogonal Hilbert transform filter banks and wavelets[END_REF]. Unfortunately, the resulting filters are not FIR and we subsequently realize causal FIR approximations of the dual filters. These approximations are obtained by minimizing the L 2 norm of the error w.r.t. the desired frequency responses. As the quality of the approximations obviously depends on the filter lengths, we introduced a 2 sample delay in the original filters so as to allow more than 16 coefficients for the dual ones. Further constraints were added in the design process, namely the desired symmetry properties and vanishing moment conditions. The latter conditions are useful in order to introduce regularity conditions on the approximate dual wavelets. Recall that a necessary and sufficient condition for the wavelets to have K vanishing moments is [START_REF] Steffen | Theory of regular M-band wavelet bases[END_REF]:

G 0 (ω) = 1 -e -ıMω 1 -e -ıω K Q(ω) , (14)
where Q(ω) is a polynomial function of e -ıω . This relation is also equivalent to the following constraints for the sequences

(g m [k]) k∈Z , m > 1: ∀m ∈ {1, . . . , M -1}, ∀k ∈ {0, . . . , K -1}, ∑ k n k g m [k] = 0.
(15) As all the considered constraints are linear, we have handled them by projection techniques. The coefficients listed in the last four columns of Table 1 have been obtained in this way, by setting K = 2, which corresponds to the same number of vanishing moments as the primal decomposition.

One of the potential problem of the proposed approach is that it does not guarantee that the designed filter bank is perfectly para-unitary. One could therefore add para-unitarity constraints in the related optimization problem. This however appeared unnecessary in the considered simulation example, as we verified that max

m,m ′ ,ℓ ∑ k g m [k]g m ′ [k -Mℓ] -δ m-m ′ δ ℓ ≤ 3.9 10 -3 . (16)
Although the upper bound is not zero as it should be ideally, its value may be considered small enough for common signal/image analysis or denoising tasks.

Associated basis functions

Given the subband filters coefficients listed in Table 1, we constructed the functions ψ m and ψ H m , for m ∈ {0, . . . , 3}, by an iterative numerical approach using the scaling equations [START_REF] Selesnick | Hilbert transform pairs of wavelet bases[END_REF]. Fig. 2 represents the scaling function and the wavelets associated to the primal 4-band filter bank. The scaling function and the wavelets associated to the dual 4-band filter bank are plotted in Fig. 3. We observe that the constructed dual wavelets possess regularity and satisfy the symmetry properties stated in Proposition 3. Symmetry axes are represented by dotted lines.

CONCLUSIONS

In this paper, we have proposed an extension of the work [START_REF] Selesnick | Hilbert transform pairs of wavelet bases[END_REF] by I. Selesnick on Hilbert transform pairs of wavelet bases to the orthogonal M-band case, by establishing phase conditions on the related filters.

We have designed filters of the dual filter bank. The resulting wavelet bases exhibit attractive symmetry and regularity properties. These encouraging design results will be adapted to forthcoming signal processing applications, especially in the field of seismic data. Taking advantages of the Hilbert pair conditions and M-band features which offer additional degrees of freedom, this new transform is promising for denoising purposes. Furthermore, the extension to 2-D or higher dimensions yields an improved directionality.

n h 0 (n) h 1 (n) h 2 (n) h 3 (n) g 0 (n) g 1 (n) g 2 (n) g 3 (n) 0 - - - - - 0

Figure 2 :

 2 Figure 2: (a) Scaling function ψ 0 and (b) wavelet ψ 1 , (c) wavelet ψ 2 , (d) wavelet ψ 3 .

Figure 3 :

 3 Figure 3: (a) Scaling function ψ H 0 and (b) wavelet ψ H 1 , (c) wavelet ψ H 2 , (d) wavelet ψ H 3 .

  Let the conditions established in Section 2.3 be satisfied. If the low-pass impulse response (h 0 [k]) k∈Z is symmetric w.r.t. (L -1)/2 where L ∈ Z and, for m ∈ {1, . . .

, M -1}, (h m [k]) k∈Z is symmetric (resp. antisymmetric) w.r.t. (L -1)/2, then (g 0 [k]) k∈Z is symmetric w.r.t. (L + M)/2 -1 and, for m ∈ {1, . . . , M -1}, (g m [k]

) k∈Z is antisymmetric (resp. symmetric) w.r.t. L/2 -1.

Table 1 :

 1 Coefficients of the synthesis filter banks for M = 4; half the coefficients are provided: h 0 , h 2 (resp. h 1 , h 3 ) are symmetric (resp. antisymmetric) of length 16; g 0 is symmetric of length 23, g 1 , g 3 are symmetric of length 19 and g 2 is antisymmetric of length 19.

	.0039 -0.0017 -0.0078 0.0043


J.-F Bercher

Tsallis distribution as a standard maximum entropy solution with 'tail' constraint

Keywords: Maximum entropy principle, Rényi entropy, nonextensive statistics PACS: 89.70.cf, 02.50.-r, 05.90.+m

We show that Tsallis' distributions can be derived from the standard (Shannon) maximum entropy setting, by incorporating a constraint on the divergence between the distribution and another distribution imagined as its tail. In this setting, we find an underlying entropy which is the Rényi entropy. Furthermore, escort distributions and generalized means appear as a direct consequence of the construction. Finally, the "maximum entropy tail distribution" is identified as a Generalized Pareto Distribution.

Introduction

The maximizers of a special entropy, the Tsallis entropy [START_REF] Tsallis | Possible generalization of Boltzmann-Gibbs statistics[END_REF], with suitable constraints, are often called Tsallis distributions. It is worth mentioning that the maximization of any monotoneous transform of Tsallis entropy, with the same constraints, leads to the same maximizers. This is in particular the case of Rényi entropy. In applied fields, Tsallis distributions (q-distributions) have encountered a large success because of their remarkable agreement with experimental data, see for instance [START_REF] Beck | Generalized statistical mechanics of cosmic rays[END_REF][START_REF] Tsallis | Entropic nonextensivity: a possible measure of complexity[END_REF] and references therein. These distributions are of very high interest in many physical systems, since they can exhibit heavy-tails, and model power-law phenomena. Indeed, power-laws are especially interesting since they appear widely in physics, biology, economy, and many other fields [START_REF] Newman | Power laws, Pareto distributions and Zipf's law[END_REF]. Tsallis distributions are similar to Generalized Pareto Distributions, which also have an high interest in other fields, namely reliability theory [START_REF] Asadi | The mean residual life function of a k-out-of-n structure at the system level[END_REF], climatology [START_REF] Montfort | Generalized Pareto distribution applied to rainfall depths[END_REF], radar imaging [START_REF] Lacour | Statistical characterization of active sonar reverberation using extreme value theory[END_REF] or actuarial sciences [START_REF] Cebrian | Generalized Pareto fit to the society of actuaries' large claims database[END_REF]. Hence, a remarkable feature of the maximum Tsallis entropy construction is its ability to exhibit heavy tailed distributions. Furthermore, the essence of Pickands' extreme values theorem [START_REF] Pickands | Statistical inference using extreme order statistics[END_REF] is that the distribution of excesses over a threshold converges, under wide conditions, to a q-exponential [START_REF] Bercher | A new look at q-exponential distributions via excess statistics[END_REF]. Following the idea that an interest of Tsallis distributions is in fact a tale of tails, in the words of [START_REF] Montroll | Maximum entropy formalism, fractals, scaling phenomena, and 1/f noise: A tale of tails[END_REF], we suggest that this kind of distributions can also be obtained from the familar (Shannon) maximum entropy setting, by the introduction of an appropriate constraint. In this setting, we show that Rényi entropy appears naturally, and also obtain a natural interpretation of nonextensive escort distributions, 'generalized means' and entropic index.

Email address: jf.bercher@esiee.fr (J.-F. Bercher). 1 On sabbatical leave from ESIEE-Paris, France.

Maximum entropy with 'tail' constraint

Jaynes' maximum entropy principle [START_REF] Jaynes | Information theory and statistical mechanics[END_REF][START_REF] Jaynes | On the rationale of maximum entropy methods[END_REF] suggests that the least biased probability distribution that describes a partially-known system is the probability distribution with maximum entropy compatible with all the available prior information. The Kullback-Leibler information divergence D(P ||Q) measures the divergence of a distribution P to another distribution Q.

Problem setting

We call 'tail distribution' the probability density function (pdf) of the excesses X u of a variate X over a threshold u. It is not only the tail part of the distribution, but is shifted by u toward the origin, and, as a pdf, is normalized to 1. Let now P and Q be two probability density functions. As a guideline, we find useful to imagine P as the tail distribution of Q: in such a case, these distributions are closely related. Our idea is to account more generally for an existing relationship between two distributions. Therefore, we propose to parametrize the relation between a candidate distribution and its 'parent' by the value θ of the Kullback-Leibler divergence between them: D(P ||Q) = θ. This define a set of possible distributions, and the solution is selected in this set, according to the maximum entropy principle, as the distribution P with maximum entropy. This writes as follows:

   max P H(P ) = -P (x) log P (x)dx s.t. D(P ||Q) = P (x) log P (x) Q(x) dx = θ. (1)
The definition of D(P ||Q) requires that P is absolutely continuous with respect to Q. It is understood, as usual, that 0 log 0 = 0 log 0/a = 0 log 0/0 = 0 and that a log(a/0) = +∞ if a > 0.

An alternative formulation, which leads to the same solution, could be to look for the distribution with minimum divergence to Q, in the set of all distributions with a given entropy:

   min P D(P ||Q) s.t. H(P ) = θ ′ .
(

) 2
It is clear that the optimum distribution arising from this procedure will not be the 'exact' tail distribution, since the distribution of the excesses X u = X -u|X > u, which reads X u = X -u if X > u (conditioned variable), has the pdf P u (x) = Q(x + u)/ Q(u), for x ≥ 0, and where Q(u) = +∞ u Q(x)dx is the so-called survival function. It is rather the maximum entropy variant, when the "tail distribution" P is constrained to be at given divergence θ to the parent distribution Q. In fact, we do not need to rely on tails in this construction, but we simply introduce a constraint on the Kullback-Leibler divergence to a 'parent' distribution in order to account for a general relationship between these distributions.

The maximum entropy solution and first consequences

The solution of (1) is easily derived using standard maximum entropy results, e.g. [START_REF] Csiszár | Why least squares and maximum entropy? an axiomatic approach to inference for linear inverse problems[END_REF]. In this case, we obtain

P (x) = 1 Z(λ) e λ log P (x) Q(x) (3)
with λ the Lagrange parameter associated to the constraint D(P ||Q) = θ, and Z(λ) the partition function. It can also be reduced to

P (x) = 1 Z(λ) 1 1-λ Q(x) λ λ-1 (4) with log Z(λ) = (1 -λ) log Q(x) λ λ-1 dx. (5)
As usual in the maximum entropy setting, the constraint and partition function are linked by

d log Z(λ) dλ = d dλ log e λ log P (x) Q(x) dx = θ. (6)
Let us now denote q = λ λ -1 .

Then, it clearly appears on the one hand that the 'maximum entropy tail' (4), which can be rewritten as

P (x) = Q(x) q Q(x) q dx , (8)
is the escort distribution of nonextensive thermostatistics. On the other hand, the log-partition function [START_REF] Asadi | The mean residual life function of a k-out-of-n structure at the system level[END_REF] becomes

log Z(q) = 1 1 -q log Q(x) q dx (9)
where the right-hand side has exactly the form of the Rényi entropy of distribution Q. It shall be mentioned that the Rényi entropy has no definite concavity for q > 1 and is not Lesche stable [START_REF] Lesche | Instabilities of Rényi entropies[END_REF][START_REF] Abe | Stability of Tsallis entropy and instabilities of rényi and normalized Tsallis entropies: A basis for q-exponential distributions[END_REF]. Its use in statistical physics has been discussed, e.g. [START_REF] Jizba | Observability of Rényi's entropy[END_REF]. In our context, it appears as a by-product of our original problem (1) which involves the maximization of the standard Boltzmann-Shannon entropy.

For the optimum distribution, we can also observe that the maximum (Shannon) entropy reduces to

H(P ) = -P (x) log P (x)dx (10) = -P (x) λ log P (x) Q(x) -log Z(λ) dx (11)
= -λ P (x) log

P (x) Q(x) dx + log Z(λ) (12)
= -λθ + log Z(λ) (13)
where the last relation is obtained using the definition of θ in [START_REF] Tsallis | Possible generalization of Boltzmann-Gibbs statistics[END_REF]. Therefore, we obtain that the maximum entropy problem (1) has for optimum value the Rényi entropy (9) with index q, minus a linear function of the constraint.

Solution with an additional observation constraint

Suppose now that the original problem is completed by an additional constraint. Indeed, the definition of equilibrium distributions usually need to take into account observation constraints. For instance, one often has to account for an observable defined as a mean value under distribution P . This is very classical in maximum entropy approaches. The constraint writes m

= E P [X] . (14)
Since P is the escort distribution (8) of Q, this mean constraint is also the generalized mean constraint of nonextensive thermostatistics:

m = xP (x)dx = 1 Q(x) q dx xQ(x) q dx. ( 15
)
When P is thought as the 'tail', m is the well-known mean residual lifetime (in reliability theory), or expected future lifetime (in survival analysis).

Because of the additional constraint, determination of the 'maximum entropy tail' distribution P amounts to further maximize the entropy in [START_REF] Jaynes | On the rationale of maximum entropy methods[END_REF] subject to that constraint. Using (9), we obtain

   max Q 1 1-q log Q(x) q dx s.t. E P [X] = m, (16)
which leads to

Q(x) ∝ (1 -β(1 -q)(x -m)) 1 1-q (17)
for the value of the parameter β such that E P [X] = m, with P given by [START_REF] Cebrian | Generalized Pareto fit to the society of actuaries' large claims database[END_REF]. As far as the latter is concerned, its expression is simply

P (x) ∝ (1 -β(1 -q)(x -m)) q 1-q . ( 18
)
This relation can also be rearranged as

P (x) ∝ (1 + β(q -1)(x -m)) -1 q-1 -1 , (19)
which is exactly in the form of the Generalized Pareto Distribution [START_REF] Johnson | Continuous Univariate Distributions[END_REF], with shape factor (q -1) and scale factor β. Using [START_REF] Lacour | Statistical characterization of active sonar reverberation using extreme value theory[END_REF], the exponent in (19) reduces to -λ, so that the distribution asymptotically behaves as a power-law with exponent -λ.

Observe that with θ = 0 in (1), we readily have P = Q and q = 1, and P in (19) reduces to the classical Boltzmann-Gibbs canonical distribution. As a final comment, let us note that the maximum entropy tail distribution is stable with respect to thresholding, in the sense that the result remains in the same family, with entropic index q but a different scale parameter. If u is the threshold, then the pdf of excesses over u, with P the parent distribution, is

P u (x) ∝ P (x + u) ∝ (1 + β(q -1)(x + u -m)) -1 q-1 -1 ∝ (1 + β ′ (q -1)(x -m)) -1 q-1 -1 , (20)
where the last term, with β ′ = β/(1 + βu), is obtained by factoring (1 + βu). This highlights the particular status of the 'maximum entropy tail distribution' as a tail distribution.

Conclusion

In this Letter, we followed two ideas. First, that an important feature of Tsallis' distributions is their ability to model the tail of distributions, particularly of those with heavy tails. The second idea is that these distributions should appear in a standard maximum entropy setting. This led us to the introduction of a constraint on the divergence between two distributions, one of them being imagined as the tail distribution. This constraint accounts for an existing general relationship between two distributions. We showed that within this construction, the escort distributions and generalized means of nonextensive statistics appear very naturally. We also obtained that the original maximum (Shannon) entropy reduces to a maximum Rényi entropy; or equivalently to the maximization of Tsallis entropy. As far as the entropic index q is concerned, it is simply associated to the value of the divergence between the distribution Q and its escort distribution P . Finally, the 'maximum entropy tail' distributions, in the sense adopted in this Letter, are found to be Generalized Pareto Distributions, which have proved very useful for modeling heavy-tailed distributions in many applied problems. We believe that this construction can be useful to workers in the field.


Nelly Pustelnik
email: nelly.pustelnik@univ-paris-est.fr

Caroline Chaux
email: caroline.chaux@univ-paris-est.fr

Jean-Christophe Pesquet
email: jean-christophe.pesquet@univ-paris-est.fr

A WAVELET-BASED QUADRATIC EXTENSION METHOD FOR IMAGE DECONVOLUTION IN THE PRESENCE OF POISSON NOISE

Keywords: Deconvolution, iterative methods, Poisson distributions, optimization methods, wavelet transforms

Iterative optimization algorithms such as the forward-backward and Douglas-Rachford algorithms have recently gained much popularity since they provide efficient solutions to a wide class of non-smooth convex minimization problems arising in signal/image recovery. However, when images are degraded by a convolution operator and a Poisson noise, a particular attention must be paid to the associated minimization problem. To solve it, we propose a new optimization method which consists of two nested iterative steps. The effectiveness of the proposed method is demonstrated via numerical comparisons.

INTRODUCTION

In many application areas, image restoration plays a primary role. When data are degraded by a convolutive blur and the addition of Gaussian noise, a large panel of methods can be used for restoration purposes. However, in real applications, the nature of the encountered noise often differs from a Gaussian one. For example, Poisson noise (which is dependent on image intensity) is found in tomography, astronomy and microscopy. This kind of noise may appear much more challenging to cope with. In this work, we consider the following degradation model:

z = P α (T y) ( 1
)
where y is the original image degraded by a blur operator T and contaminated by a Poisson noise with scaling factor α.

The effect of the noise is denoted by P α .

The simplest way to take into account Poisson noise is to process the data as if they were corrupted with Gaussian noise and consequently, to apply standard restoration tools [START_REF] Bouman | A unified approach to statistical tomography using coordinate descent optimization[END_REF][START_REF] Verhaeghe | Dynamic PET reconstruction using wavelet regularization with adapted basis functions[END_REF]. A better approach consists of applying a pre-processing (also called Variance Stabilizing Transform) on the data. For example, the Anscombe transform [START_REF] Anscombe | The transformation, of Poisson, binomial and negative-binomial data[END_REF] can be used to pre-process Poisson data and then, to restore the transformed data as if they were corrupted by Gaussian noise. Fryzlewicz and Nason explored another stabilizing variance approach in [START_REF] Fryzlewicz | A Haar-Fisz algorithm for poisson intensity estimation[END_REF] named the Haar-Fisz transform. The main drawback of this pre-processing lies in the fact that it uses the Haar decomposition which is known to be especially efficient for blockwise constant images. More recently, iterative algorithms based on a variational approach were proposed in the Poisson denoising context. In [START_REF] Sardy | Automatic smoothing with wavelets for a wide class of distributions[END_REF], Sardy et al. developed an ℓ 1 -penalized likelihood approach and more recently, in [START_REF] Combettes | A Douglas-Rachford splitting approach to nonsmooth convex variational signal recovery[END_REF], the authors used the Douglas-Rachford algorithm within a convex optimization framework. For deconvolution problems, algorithms based on the minimization of the Kullback-Leibler divergence were studied in [START_REF] Byrne | Iterative image reconstruction algorithms based on crossentropy minimization[END_REF]. An extension of this approach was investigated in [START_REF] Dey | Richardson-Lucy algorithm with total variation regularization for 3D confocal microscope deconvolution[END_REF] to incorporate a total variation penality term. More recently, a block iterative method for binary tomographic reconstruction based on convex contraint sets was proposed in [START_REF] Capricelli | A convex programming algorithm for noisy discrete tomoraphy[END_REF]. In [START_REF] Dupé | Image deconvolution under Poisson noise using sparse representations and proximal thresholding iteration[END_REF][START_REF] Pustelnik | A constrained forwardbackward algorithm for image recovery problems[END_REF], the authors considered algorithms mixing forward-backward [START_REF] Chaux | A variational formulation for frame-based inverse problems[END_REF] and Douglas-Rachford steps. On the one hand, Dupé et al. [START_REF] Dupé | Image deconvolution under Poisson noise using sparse representations and proximal thresholding iteration[END_REF] investigated an elegant adaptation of the Anscombe approach. On the other hand, a quadratic extension dealing with signal dependent Gaussian noise was introduced in [START_REF] Pustelnik | A constrained forwardbackward algorithm for image recovery problems[END_REF]. Note that a different polynomial approximation of the objective function was proposed in [START_REF] Fessler | Hybrid poisson/polynomial objective functions for tomographic image reconstruction from transmission scans[END_REF]. The objective of this paper is to adapt the extension proposed in [START_REF] Pustelnik | A constrained forwardbackward algorithm for image recovery problems[END_REF] in order to deal with Poisson distributions and, to investigate its properties. Notice that the proposed method relies on a wavelet-like representation as in some of the aforementioned approaches.

The paper is organized as follows: in Section 2, we formulate the restoration problem from a frame representation perspective. A Bayesian interpretation of the related convex minimization problem is given. Then, we propose a quadratic extension technique to circumvent the problem of the non Lipschitz differentiability of the likelihood of the Poisson data. In Section 3, we present the new algorithm combining a forward-backward and a Dykstra's step so as to restore blurred images contaminated by Poisson noise. Finally, the effectiveness of the proposed approach is demonstrated via simulation results.

IMAGE DECONVOLUTION IN THE PRESENCE OF POISSON NOISE

Problem statement

We consider the degradation model ( 1). An efficient probabilistic prior model on the unknown image y is adopted by considering a frame representation of this image [START_REF] Chaux | A variational formulation for frame-based inverse problems[END_REF]. We thus use a linear representation of the form y = F * x where 1) can be reformulated as z = P α (T F * x) where x is the unknown vector of frame coefficients.

F * : H → G is a tight frame synthesis operator with H = R K , G = R N and K ≥ N . This means that F * • F = νId with ν ∈]0, +∞[. In this context, (

Bayesian framework

The image u = T y (resp. z) is viewed as a realization of a nonnegative real-valued random vector U = (U

(i) ) 1≤i≤N (resp. Z = (Z (i) ) 1≤i≤N ). Conditionally to U = (u (i) ) 1≤i≤N
∈ G, the random vector Z is assumed to have independent components with conditional probability: ∀z (i) ∈ N,

P(Z (i) = z (i) | U (i) = u (i) ) = (αu (i) ) z (i) z (i) ! exp(-αu (i) ) (2)
where α ∈]0, +∞[ is a scaling parameter.

We then assume that the vector x of frame coefficients is a realization of a random vector X with independent components and, each component

X (k)
1≤k≤K has a probability density given by: ∀x

(k) ∈ R, p X (k) (x (k) ) ∝ exp(-φ k (x (k) ))
where φ k is a finite convex function allowing us to model frame coefficients. A large choice of potential functions (φ k ) 1≤k≤K is possible as shown in [START_REF] Chaux | A variational formulation for frame-based inverse problems[END_REF]. Under these assumptions, it is readily shown that a MAP estimator of the vector of frame coefficients x can be obtained from z by solving:

min x=(x (k) ) 1≤k≤K ∈H K X k=1 φ k (x (k) ) -ln P(Z = z | U = T F * x).

Link with convex optimization

One generic problem in convex optimization is to determine:

min x∈C f (x) + g(x) (3)
where C is a closed convex subset of H, f and g are functions in Γ 0 (H). Γ 0 (H) denotes the class of lower semicontinuous convex functions taking their values in ]-∞, +∞].

In our restoration problem, g denotes the fidelity term related to the Poisson distribution of the noise and f the a priori term on x. Thus, g is chosen as, ∀x ∈ H, g(x) = Ψ(T F * x), where

∀u = u (i) 1≤i≤N ∈ G Ψ(u) = N i=1 ψ i u (i) . (4)
From ( 2), it can be deduced that: ∀υ ∈ R,

ψ i (υ) =        αυ -z (i) + z (i) ln z (i) αυ if z (i) > 0 and υ > 0, αυ if z (i) = 0 and υ ≥ 0, +∞ otherwise.
(5) Furthermore, the prior f is defined as:

(∀x = x (k) 1≤k≤K ∈ R K ) f (x) = K k=1 φ k x (k) . (6)
Finally, as it is customary in image processing, a constraint on the range of the pixel values is introduced by setting

C = x ∈ H F * x ∈ [0, 255] N . (7)

Lipschitz differentiability issue

The functions f and g as defined above are in Γ 0 (H) but an additional assumption is often required to ensure the convergence of optimization algorithms such as the modified forward-backward algorithm in [START_REF] Pustelnik | A constrained forwardbackward algorithm for image recovery problems[END_REF], namely g should be β-Lipschitz differentiable on C with β ∈]0, +∞[. Unfortunately the assumption is not satisfied by g.

In [START_REF] Pustelnik | A constrained forwardbackward algorithm for image recovery problems[END_REF], a quadratic extension was proposed so as to improve convergence profiles. In the present case, we will show that this technique can be adapted in order to obtain a lower approximation g θ of g, which is Lipschitz differentiable on C. This function is defined as:

g θ = Ψ θ • T • F * , where θ ∈ ]0, +∞[ and Ψ θ is given by ∀u = u (i) 1≤i≤N ∈ G Ψ θ (u) = N i=1 ψ θ,i u (i) . ( 8
)
The functions (ψ θ,i ) 1≤i≤N are here defined as: ∀υ ∈ R,

ψ θ,i (υ) =        θ 2 υ 2 + ζ i,1 (θ) υ + ζ i,0 (θ) if z (i) > 0 and -ǫ(θ) ≤ υ < υi(θ) αυ if z (i) = 0 and -ǫ(θ) ≤ υ < 0 ψ i (υ) otherwise,
where

υ i (θ) = (z (i) /θ) 1/2 is such that 0 ≤ ψ ′′ i (υ) ≤ θ ⇔ υ ≥ υ i (θ),
and the constants ζ i,0 (θ) and ζ i,1 (θ) are chosen so as to ensure the continuity of ψ θ,i in υ i (θ). The validity of this approximation is secured by the following result: Proposition 2.1 [START_REF] Chaux | Nested iterative algorithms for convex constrained image recovery problems[END_REF] Assume that (i) ǫ is a decreasing positive function such that lim θ→+∞ ǫ(θ) = 0,

(ii) T [0, 255] N ⊂ [0, +∞[ N , (iii) f is coercive or C is bounded, (iv) f is strictly convex on C.
Then, there exists θ ∈ ]0, +∞[ such that, for every θ ∈ [θ, +∞[, f + g θ has a unique minimizer on the convex set C which is the minimizer of f + g on C.

The proposed quadratic extension is illustrated in Fig. 1. When θ becomes greater than some lower value θ, the solution of the original minimization problem is perfectly found as g is approximated by g θ . The choice of θ will be discussed in Section 4.

PROPOSED ALGORITHM

Optimization background

The approximated version of Problem (3) can be expressed as:

min x∈H f (x) + g θ (x) + ι C (x) (9)
where ι C denotes the indicator function of C, i.e., ∀x ∈ H, ι C (x) = 0 if x ∈ C, +∞ otherwise. Among convex optimization methods, the forward-backward algorithm appears as an appropriate solution for minimizing (9) when f is non necessarily smooth. It is based on a useful tool called the proximity operator. We recall that the proximity operator of

f ∈ Γ 0 (H) is prox f : H → H : x → prox f x where prox f x is the unique minimizer of f + 1 2 . -x 2 .
This operator generalizes the notion of projection as prox ιC = P C where P C denotes the projection onto C.

Forward-backward approach

Let x 0 ∈ C be an initial image. The algorithm constructs a sequence (x n ) n≥1 by the iteration: for every n ∈ N,

x n+1 = x n + λ n prox ιC +γnf (x n -γ n ∇g θ (x n )) -x n
where γ n ∈]0, 2(νθ T 2 ) -1 [ is the algorithm step-size and λ n ∈]0, 1[ is a relaxation parameter. Under these conditions, it is known that, if (∀n ∈ N) x n ∈ C, then (x n ) n∈N converges to a solution to Problem [START_REF] Capricelli | A convex programming algorithm for noisy discrete tomoraphy[END_REF] [START_REF] Chaux | A variational formulation for frame-based inverse problems[END_REF]Theorem 5.4].

Dykstra's algorithm

A main difficulty in the application of the forward-backward algorithm is the determination of prox ιC +γnf [START_REF] Pustelnik | A constrained forwardbackward algorithm for image recovery problems[END_REF]. To perform this task, we propose to use Dykstra's algorithm [START_REF] Bauschke | A Dykstra-like algorithm for two monotone operators[END_REF]. Set r 0 ∈ H and p 0 = q 0 = 0. Sequences (r m ) m≥1 and (s m ) m≥1 are generated by the following routine: ∀m ∈ N,

(s m , p m+1 , r m+1 , q m+1 ) = D f (r m , p m , q m ) ⇔          s m = P C (r m + p m ) p m+1 = r m + p m -s m r m+1 = prox f (s m + q m ) q m+1 = s m + q m -r m+1 .
Then (r m ) m∈N and (s m ) m∈N converge to the solution of prox f +ιC (r 0 ) [15, Theorem 3.3].

Resulting algorithm

The proposed method consists of the following combination of the two previous algorithms:

Set x 0 ∈ C and n = 0. Set r n,0 = x n -γ n ∇g θ (x n ) and p n,0 = q n,0 = 0. For m = 0, . . . , M n a) (s n,m , p n,m+1 , r n,m+1 , q n,m+1 ) = D γnf (r n,m , p n,m , q n,m ) b) If s n,m = s n,m-1 , goto . Set x n+1 = x n + λ n s n,m -x n .
Increment n (n ← n + 1) and goto .

The following result can then be deduced from the convergence properties of the forward-backward and Dykstra's algorithms.

Proposition 3.1 Under Assumptions (ii) and (iii) of Proposition 2.1, there exists a sequence of positive integers (M n ) n∈N such that, if (∀n ∈ N) M n ≥ M n then, (x n ) n∈N converges to a solution to Problem (9).

Remark

Other convex optimization algorithms can be envisaged to solve Problem [START_REF] Capricelli | A convex programming algorithm for noisy discrete tomoraphy[END_REF]. On the one hand, it is possible to replace the Dykstra's inner loop by iterations of the Douglas-Rachford algorithm [START_REF] Combettes | A Douglas-Rachford splitting approach to nonsmooth convex variational signal recovery[END_REF]. On the other hand, a minimization strategy relying on the Douglas-Rachford algorithm can be employed where an inner loop consists of iterations of the forward-backward algorithm [START_REF] Chaux | A variational formulation for frame-based inverse problems[END_REF]. In our simulations however, it appeared that the algorithm we propose in this paper is more effective in terms of convergence speed than these two alternatives.

SIMULATION RESULTS

The objective of the presented experiment is to recover a microscopy image from a degraded observation z generated according to Model (1). The degradation is obtained by applying a 5 × 5 positive uniform blur with T = 1. To restore this image, a tight frame version of the dual-tree transform (DTT) proposed in [START_REF] Chaux | Image analysis using a dualtree M -band wavelet transform[END_REF] In Table 1, we provide a quantitative comparison between the regularized Expectation-Maximization algorithm proposed in [START_REF] Byrne | Iterative image reconstruction algorithms based on crossentropy minimization[END_REF], the Haar-Fisz transform [START_REF] Fryzlewicz | A Haar-Fisz algorithm for poisson intensity estimation[END_REF], the Anscombe approach in [START_REF] Dupé | Image deconvolution under Poisson noise using sparse representations and proximal thresholding iteration[END_REF] and the proposed method for different values of the quadratic extension parameter θ. Table 1. SNR for the microscopy image. For quadratic extension method, iteration number <500 when θ = 0.001 and ≃ 1000 when θ = 0.1.

As mentioned in Section 2.4, a large value of θ allows us to improve the restoration performance (closer approximation g θ of g), whereas a too small value may yield poorly restored images. However, the step-size γ n being inversely proportional to θ, a compromise has to be done between the convergence rate of the algorithm and the accuracy to the model. It can be noticed that a significant gain is obtained in comparison with other methods except for the case when α = 1 for which the Anscombe approach gives similar results.

CONCLUSIONS

A new restoration algorithm has been proposed for waveletbased restoration in the presence of Poisson noise. This algo-rithm allows us to minimize a wide class of convex objective functions where the fidelity term is directly related to the Poisson likelihood. A quadratic extension is necessary to ensure the gradient Lipschitz property of the smooth term. Although not thoroughly discussed in this paper, this extension is also beneficial to the improvement of the convergence rate.

Fig. 1 .

 1 Fig. 1. Quadratic approximation: Poisson likelihood (continuous line) and its extensions for θ = 0.3 (dashed line), θ = 0.07 (dash-dot line) and θ = 0.01 (dotted line).

  (ν = 2) using Alkin and Caglar 4-band filter bank has been employed. Strictly convex non-differentiable potential functions φ k are chosen, of the form ω k | . | p k + χ k | . | where (ω k , χ k ) ∈ ]0, +∞[ 2 and p k ∈ {4/3, 3/2, 2}. We have used the algorithm proposed in Section 3.4 with x 0 = P C z (see (7) for the definition of C) and γ n = 0.995 θ -1 . Visual results are shown in the following figure. Original Degraded, α = 0.1 Restored, θ = 0.1


Nelly Pustelnik
email: nelly.pustelnik@univ-paris-est.fr

Caroline Chaux
email: caroline.chaux@univ-paris-est.fr

Jean-Christophe Pesquet
email: jean-christophe.pesquet@univ-paris-est.fr

HYBRID REGULARIZATION FOR DATA RESTORATION IN THE PRESENCE OF POISSON NOISE

During the last five years, several convex optimization algorithms have been proposed for solving inverse problems. Most of the time, they allow us to minimize a criterion composed of two terms one of which permits to "stabilize" the solution. Different choices are possible for the so-called regularization term, which plays a prominent role for solving ill-posed problems. While a total variation regularization introduces staircase effects, a wavelet regularization may bring other kinds of visual artefacts. A compromise can be envisaged combining these regularization functions. In the context of Poisson data, we propose in this paper an algorithm to achieve the minimization of the associated (possibly constrained) convex optimization problem.

INTRODUCTION

Different algorithms in convex optimization have been proposed to minimize regularized criteria used to provide solutions to ill-posed problems. In this paper, we are more specifically interested in restoring degraded images when the noise is Poisson distributed.

A way to take into account Poisson noise consists of applying a pre-processing on the data (like the Anscombe [START_REF] Anscombe | The transformation, of Poisson, binomial and negative-binomial data[END_REF] or the Haar-Fisz [START_REF] Fryzlewicz | A Haar-Fisz algorithm for poisson intensity estimation[END_REF] transforms) in order to stabilize the variance and then to apply standard restoration tools [START_REF] Bouman | A unified approach to statistical tomography using coordinate descent optimization[END_REF][START_REF] Verhaeghe | Dynamic PET reconstruction using wavelet regularization with adapted basis functions[END_REF].

One of the first restoration method aiming at maximizing the Poisson likelihood is the Richardson-Lucy algorithm [START_REF] Richardson | Bayesian-based iterative method of image restoration[END_REF][START_REF] Lucy | An iterative technique for the rectification of observed distributions[END_REF]. However, the main drawback of this method is that it does not allow us to incorporate information about the target solution, except the positivity constraint in an implicit manner. An improved form of this algorithm involving a total variation [START_REF] Rudin | Nonlinear total variation based noise removal algorithms[END_REF] regularization was recently proposed in [START_REF] Dey | Richardson-Lucy algorithm with total variation regularization for 3D confocal microscope deconvolution[END_REF]. Due to the equivalence between the maximization of the Poisson likelihood and the mimization of a Kullback-Leibler divergence, other forms of multiplicative iterative algorithms have been proposed in order to take into account specific forms of regularizations [START_REF] Byrne | Iterative image reconstruction algorithms based on crossentropy minimization[END_REF].

However, during the last decade, much interest has been gained in introducing a priori information about the target image in a transformed domain. Actually, wavelets are often used due to their ability to provide sparse representations for many classes of regular signals. In this respect, redundant frames (overcomplete transforms) constitute more flexible tools than orthonormal bases for building linear representations of images. A number of recent works [START_REF] Chaux | A variational formulation for frame-based inverse problems[END_REF] have emphasized the interest in using specific tight frames for performing geometrical analyses of images.

Existing works for solving inverse problems in a frame are often restricted to noise with log-likelihood having a Lip-schitz continuous gradient [START_REF] Chaux | A variational formulation for frame-based inverse problems[END_REF]. To alleviate this restriction, the authors in [START_REF] Dupé | Image deconvolution under Poisson noise using sparse representations and proximal thresholding iteration[END_REF][START_REF] Chaux | Nested iterative algorithms for convex constrained image recovery problems[END_REF] considered algorithms mixing forward-backward [START_REF] Chaux | A variational formulation for frame-based inverse problems[END_REF] and Douglas-Rachford algorithms. First, Dupé et al. [START_REF] Dupé | Image deconvolution under Poisson noise using sparse representations and proximal thresholding iteration[END_REF] investigated an elegant adaptation of the Anscombe approach and then, in [START_REF] Chaux | Nested iterative algorithms for convex constrained image recovery problems[END_REF] a quadratic extension technique was introduced.

One of the drawbacks of the approaches based on wavelet representations is that they may introduce visual artefacts, e.g. some lack of regularity in homogeneous areas or ringing artefacts along edges. Alternative solutions based on the use of the total variation can be employed but they often lead to so-called staircase effects. To circumvent these problems and take advantage of both approaches, authors have suggested in [START_REF] Combettes | A proximal decomposition method for solving convex variational inverse problems[END_REF][START_REF] Bioucas-Dias | An iterative algorithm for linear inverse problems with compound regularizers[END_REF] to jointly take into account the ℓ 1 -norm of the wavelet coefficients and the total variation penalization, which resulted in a significant improvement of the visual restoration quality. These works are however restricted to data corrupted by Gaussian noise.

The purpose of this paper is to extend the work in [START_REF] Combettes | A proximal decomposition method for solving convex variational inverse problems[END_REF] so as to solve linear inverse problems in the presence of Poisson noise. Let H be a real separable Hilbert space corresponding to the frame coefficient space. We aim at finding

min x∈C g(x) + f (x) + h(x) (1)
where g, f and h are assumed to be in the class Γ 0 (H ) of lower semicontinuous convex functions on H taking their values in ] -∞, +∞] which are proper (i.e. not identically equal to +∞). h is related to a total variation penalization and C is a nonempty closed convex subset of H . From a Bayesian interpretation when C = H , f + h can be viewed as an a priori term on the frame coefficients of the original image, g corresponds to the Poisson anti log-likelihood (fidelity term) and the main difficulty is to minimize (1) when the degradation is a combined effect of a linear operator (e.g. a blur) and a Poisson noise. These three functions will be detailled in Section 2.

The proposed method is based on recent developments in convex optimization. Its main advantage is that it does not require any approximation of the noise distribution unlike the Anscombe approach [START_REF] Dupé | Image deconvolution under Poisson noise using sparse representations and proximal thresholding iteration[END_REF] or the quadratic extension in [START_REF] Chaux | Nested iterative algorithms for convex constrained image recovery problems[END_REF]. The main contribution of our work consists of taking advantage of the structural properties of the linear degradation operator so as to design an efficient algorithm.

In Section 2, we formulate the minimization problem associated with the restoration of data corrupted by Poisson noise. A hybrid regularization combining a wavelet sparsity promoting term and a total variation term is considered.

The optimization algorithm used to solve this problem is also described based on the work in [START_REF] Combettes | A proximal decomposition method for solving convex variational inverse problems[END_REF]. In Section 3, we briefly recall the definition of the proximity operator which is the building block of this algorithm. The main difficulty to be addressed here is that the proximity operator of the Kullback-Leibler divergence term does not have an explicit expression. A way of solving this problem is then proposed. Finally, in Section 4, we provide numerical examples as well as algorithm implementation details.

PROBLEM FORMULATION

The degradation model is the following,

z = P α (T ȳ) (2)
where ȳ ∈ R N + is the original data of size N degraded by a non-negative valued convolutive operator T and contaminated by a Poisson noise with scaling parameter α ∈ R * + . The vector z ∈ N N ⊂ G = R N represents the observed data.

Frame representation

In inverse problems, certain physical properties of the target solution y are most suitably expressed in terms of the coefficients x = (ξ k ) k∈K⊂N (where K = {1, . . . , K}) of its representation y = ∑ k∈K ξ k e k with respect to a family of vectors (e k ) k∈K⊂N in a Hilbert space G . Recall that a family of vectors (e k ) k∈K in G constitutes a frame if there exist two constants ν and ν in ]0, +∞[ such that

(∀y ∈ G ) ν y 2 ≤ ∑ k∈K | y, e k | 2 ≤ ν y 2 .
(

The associated frame operator is the injective bounded linear operator

F : G → ℓ 2 (K) : y → ( y, e k ) k∈K , (4)
the adjoint of which is the surjective bounded linear operator

F * : ℓ 2 (K) → G : (ξ k ) k∈K → ∑ k∈K ξ k e k . (5)
When ν = ν = ν in (3), (e k ) k∈K is said to be a tight frame.

In this case, we have

F * • F = νId, (6)
where Id is the identity operator. A simple example of a tight frame is the union of m orthonormal bases, in which case ν = ν = m. Considering a frame representation, Model (2) can be re-expressed as,

z = P α (T F * x) (7)
where x represents frame coefficients of the original data (y = F * x).

Convex optimization

In the context of inverse problems, recent studies proposed to restore the original signal by solving a convex optimization problem of the form:

Find min x∈H J ∑ j=1 f j (x) (8)
where ( f j ) 1≤ j≤J are functions of Γ 0 (H ) where H = R K . In restoration methods, J is usually equal to 2. One function is a smooth term related to the observation model and the second one is a regularization term. In the case of data degraded by a Poisson noise, a standard choice for the first function is a Kullback-Leibler divergence term [START_REF] Titterington | On the iterative image space reconstruction algorithm for ECT[END_REF]. More precisely, we employ the Kullback-Leibler divergence D KL from T F * x to z, for every x ∈ H :

D KL (T F * x, z) = Ψ(T F * x) (9)
where,

∀u = (u (n) ) 1≤n≤N ∈ G , Ψ(u) = N ∑ n=1 ψ n (u (n) )
. (10) and

ψ n (υ) =        αυ -z (n) + z (n) ln z (n) αυ if z (n) > 0 and υ > 0, αυ if z (n) = 0 and υ ≥ 0, +∞ otherwise.
Concerning the regularization function, in [START_REF] Figueiredo | An EM algorithm for waveletbased image restoration[END_REF][START_REF] Daubechies | An iterative thresholding algorithm for linear inverse problems with a sparsity constraint[END_REF][START_REF] Chaux | A variational formulation for frame-based inverse problems[END_REF]] the authors consider a penalization in the wavelet domain corresponding to power functions of the coefficients with exponents p k ∈ {1, 4/3, 3/2, 2, 3, 4}. Another type of regularization that can be envisaged employs a total variation measure. Recently, in [START_REF] Combettes | A proximal decomposition method for solving convex variational inverse problems[END_REF], a hybrid regularization was proposed which gives good results in the context of an additive white Gaussian noise. We thus propose to use a similar compound regularization in the case when the data are corrupted by Poisson noise. In this context, Problem (1) is solved by setting g = Ψ • T • F * where Ψ is defined by [START_REF] Chaux | A variational formulation for frame-based inverse problems[END_REF]. The function f corresponds to the regularization term operating in the wavelet domain, which is chosen such that

∀x = (x (k) ) 1≤k≤K ∈ H , f (x) = ∑ K k=1 φ k (x (k) )
where, for ev- ery k ∈ {1, . . . , K}, φ k is a finite function of Γ 0 (R) such that lim |x (k) |→+∞ φ k (x (k) ) = +∞ and h represents a total variation term such that h = tv • F * , tv being the total variation operator defined in [13, Experiment 2]. Hence, the considered minimization problem in (1) becomes, (11) where ι C is the indicator function of a closed convex set C (for example related to support or value range constraints) such that,

min x∈H D KL (T F * x, z) + ϑ f (x) + κ tv(F * x) + ι C (x)
(∀x ∈ H ) ι C (x) = 0, if x ∈ C; +∞, otherwise. ( 12
)
Throughout the paper, it is assumed that

(TC * )∩]0, +∞[ N = ∅ with C * = F * C = {F * x | x ∈ C}.
The nonnegative real parameters ϑ and κ control the degree of smoothness in the wavelet and in the space domain, respectively. Notice that, when ϑ = 0, the regularization reduces to the standard total variation penalization and that the approach we propose in this paper thus provides also an efficient numerical solution to such a more classical penalized optimization problem.

In the class of convex optimization methods, to the best of our knowledge, only one algorithm allows us to efficiently minimize the sum of four terms, which are not differentiable. The algorithm was proposed by Combettes and Pesquet in [START_REF] Combettes | A proximal decomposition method for solving convex variational inverse problems[END_REF] and is summarized next. Algorithm 2.1 Let γ ∈ ]0, +∞[. For every j ∈ {1, . . . , J}, set (ω j ) 1≤ j≤J ∈]0, 1] J such that ∑ J j=1 ω j = 1, (y j,0 ) 1≤ j≤J ∈ H J and x 0 = ∑ J j=1 ω j y j,0 . Let (a j,l ) l∈N be a sequence in H which corresponds to possible errors in the computation of proximity operators. Then, the sequence (x l ) l≥1 is generated by the following routine: for every l ∈ N,

                    
For j = 1, . . . , J ⌊p j,l = prox γ/ω j f j y j,l + a j,l p l = ∑ J j=1 ω j p j,l λ l ∈ ]0, 2[ For j = 1, . . . , J ⌊y j,l+1 = y j,l + λ l (2 p l -x l -p j,l )

x l+1 = x l + λ l (p l -x l ) (13)
The prox operator introduced in the first loop will be defined Section 3. The sequence (x l ) l≥1 converges to a solution to Problem (8) under the following assumption.

Assumption 2.1 (i) lim x →+∞ f 1 (x) + . . . + f J (x) = +∞. (ii) ∩ J j=1 rintdom f j = ∅ (when H is finite-dimensional) 1 . (iii) ∑ l∈N λ l (2 -λ l ) = +∞. (iv) (∀ j ∈ {1, . . . , J}) ∑ l∈N λ l a j,l < +∞.
The main difficulty in applying this algorithm to our restoration problem is that it requires to compute the proximity operators associated to each of the four terms in [START_REF] Dupé | Image deconvolution under Poisson noise using sparse representations and proximal thresholding iteration[END_REF]. Closed forms of these quantities are known for the indicator function, the function f [START_REF] Chaux | A variational formulation for frame-based inverse problems[END_REF] and the total variation term (through the decomposition method proposed in [START_REF] Combettes | A proximal decomposition method for solving convex variational inverse problems[END_REF]). The main problem remains in the computation of the proximity operator of the Kullback-Leibler divergence term. In the next section, we will recall the notion of proximity operator, before providing an answer to this question.

CONVEX OPTIMIZATION TOOLS: PROXIMITY OPERATOR

Definition and properties

A fundamental tool in the study of the convex optimization methods is the proximity operator introduced by Moreau in 1962 [START_REF] Moreau | Proximité et dualité dans un espace hilbertien[END_REF]. The proximity operator of ϕ ∈ Γ 0 (H ) is defined by

prox ϕ : H → H ; x → arg min y∈H 1 2 y -x 2 + ϕ(y). ( 14
)
We thus see that prox ι C reduces to the projection P C onto the convex set C. Other examples of proximity operators corresponding to the potential functions of standard log-concave univariate probability densities have been listed in [START_REF] Chaux | A variational formulation for frame-based inverse problems[END_REF].

We now recall the proximity operator of the potential associated with a Gamma distribution which is closely related to the Kullback-Leibler divergence [START_REF] Combettes | A Douglas-Rachford splitting approach to nonsmooth convex variational signal recovery[END_REF]. 1 The relative interior of a set S of H is designated by rint S and the domain of a function f :

H →] -∞;+∞] is dom f = {x ∈ H | f (x) < +∞}.
Example 3.1 [START_REF] Chaux | A variational formulation for frame-based inverse problems[END_REF] Let ω ∈ ]0, +∞[, χ ∈ [0, +∞[, and set

ϕ : R → ]-∞, +∞] : η → -χ ln(η) + ωη, if η > 0; +∞, if η ≤ 0. ( 15
) Then, for every η ∈ R,

prox ϕ η = η -ω + |η -ω| 2 + 4χ 2 . ( 16
)
Our minimization problem being formulated in a frame representation, we will also need a property concerning the calculation of the proximity operator of the composition of a function of Γ 0 (G ) and a linear operator: Proposition 3.2 [START_REF] Combettes | A Douglas-Rachford splitting approach to nonsmooth convex variational signal recovery[END_REF] Let H and G be real Hilbert spaces, let ϕ ∈ Γ 0 (G ), and let L : H → G be a bounded linear operator. Suppose that the composition of L and L * satisfies L • L * = ς Id, for some ς ∈ ]0, +∞[. Then ϕ • L ∈ Γ 0 (H ) and

prox ϕ•L = Id + ς -1 L * • (prox ς ϕ -Id) • L. ( 17
)

Case of a convolutive operator T

In order to solve Problem ( 11), we will be interested in determining the proximity operator of g = Ψ • T • F * . As will be shown next, the proximity operator of this function can be determined in a closed form for specific cases only. However, g can be decomposed as a sum of functions for which the proximity operators can be calculated explicitly. We subsequently assume that:

Assumption 3.3 (e k ) 1≤k≤K is a tight frame of G = R N with frame constant ν > 0.
Let us now focus on function g. Let (I i ) 1≤i≤I be a partition of {1, . . . , N} in nonempty sets. For every i ∈ {1, . . . , I}, let M i be the number of elements in I i and let ϒ i : R M i → ]0, +∞[ : (η n ) n∈I i → ∑ n∈I i ψ n (η n ). We have then

g = ∑ I i=1 ϒ i • T i • F *
where T i is the linear operator from R N to R M i associated with the matrix

[t m 1 . . . t m M i ] ⊤ with I i = {m 1 , . . . , m M i }.
The following assumption will play a prominent role in the rest of paper: Assumption 3.4 For all i ∈ {1, . . . , I}, (t n ) n∈I i is a family of orthogonal vectors having the same norm σ

1/2 i with σ i > 0.
Note that this assumption is obviously satisfied when I = N, that is when ∀i ∈ {1, . . ., I}, I i reduces to a singleton.

Proposition 3.5 Under Assumptions 3.3 and 3.4, we have,

(∀i ∈ {1, . . . , I}) (∀x ∈ H ) prox ϒ i •T i •F * (x) = x + 1 νσ i FT * i (π (n) i -η (n) i ) n∈I i ( 18
)
where (η

(n) i ) n∈I i = T i F * x and (∀n ∈ I i ) π (n) i = η (n) i -ανσ i + |η (n) i -ανσ i | 2 + 4νσ i z (n) 2 . (19)

Remark

(i) This result is a consequence of Proposition 3.2 by setting L = T i F * and by using Example 3.1 to derive prox νσ i ϒ i . (ii) It can be noticed that the application of T i or T * i reduces to standard operations in signal processing. The application of T i consists of two steps: a convolution with the degradation filter and a decimation for selected locations (n ∈ I i ). The application of T * i also consists of two steps: an interpolation step (by setting zero everywhere except for indices n ∈ I i ) followed by a convolution with the filter with conjugate frequency response.

To reduce the number of proximity operators to be computed, we want to find the smallest integer I such that, for every i ∈ {1, . . . , I}, (t n ) n∈I i is an orthogonal basis. Consider the particular case of a deconvolution problem for images of size

N 1 × N 2 (N = N 1 N 2 ) involving a periodic convolution where Q 1 × Q 2 is
the kernel size. The operator T is associated with a matrix having a circulant-block circulant structure [START_REF] Golub | Matrix computations[END_REF] and consequently,

   t ⊤ 1 . . . t ⊤ N    =                  Θ 0 O ... O Θ Q 1 -1 ... Θ 1 Θ 1 Θ 0 . . . . . . . . . . . . . . . . . . . . . . . . Θ Q 1 -1 Θ Q 1 -1 . . . . . . O O . . . . . . . . . . . . . . . . . . . . . . . . O O ... O Θ Q 1 -1 ... Θ 1 Θ 0                  ( 20
)
where O is the null matrix of size N 2 × N 2 and (∀q 1 ∈ {0, . . . , Q 1 -1}) Θ q 1 is the matrix of size N 2 × N 2 defined by

Θ q 1 =                  θ q 1 ,0 0 ... 0 θ q 1 ,Q 2 -1 ... Θ q 1 ,1 θ q 1 ,1 θ q 1 ,0 0
0 θ q 1 ,Q 2 -1 ... θ q 1 ,1 θ q 1 ,0                  . (21)
(θ q 1 ,q 2 ) 0≤q 1 <Q 1 ,0≤q 2 <Q 2 denotes here the point spread function of the degradation filter.

Let us define

I 1 = min{i 1 ≥ Q 1 | N 1 = 0 mod i 1 } and, I 2 = min{i 2 ≥ Q 2 | N 2 = 0 mod i 2 }.
In order to satisfy Assumption 3.4, we subsequently set

I = I 1 I 2 .
For all (i 1 , i 2 ) ∈ {1, . . . , I 1 } × {1, . . . , I 2 }, set i = i 2 + I 2 (i 1 -1) and define

I i = {n 2 + N 2 (n 1 -1) | (n 1 , n 2 ) ∈ {1, . . . , N 1 } × {1, . . ., N 2 }
and

n 1 = i 1 mod I 1 , n 2 = i 2 mod I 2 }. (22)
Then, Assumption 3.4 holds and, for all i ∈ {1, . . . , I},

σ i = ∑ Q 1 -1 q 1 =0 ∑ Q 2 -1 q 2 =0 |θ q 1 ,q 2 | 2 .
We have therefore defined sets (I i ) 1≤i≤I for which we can compute the associated proximity operators as expressed by Proposition 3.5. Note that calculations for other forms of convolution operations (e.g. including zero-padding) can be achieved in a similar way.

Regarding Algorithm 2.1, the J proximity operators to be computed are: the I proximity operators corresponding to prox γ ω i ϒ i •T i •F * , the four proximity operators associated with the total variation term [START_REF] Combettes | A proximal decomposition method for solving convex variational inverse problems[END_REF], the proximity operator of the regularization function f and also, the projection onto the convex set C. To sum up, J = I + 6 proximity operator computations are needed in order to solve Problem [START_REF] Dupé | Image deconvolution under Poisson noise using sparse representations and proximal thresholding iteration[END_REF].

EXPERIMENTAL RESULTS

In our simulations, our objectives are twofold: we will first be interested in studying the influence of the combination of total variation and wavelet regularization terms and then, we will compare the results obtained by using the proposed algorithm with those corresponding to state-of-the-art methods. Two test images (N 1 = N 2 = 256) will be considered (see Fig. 1). For both examples, T is a uniform blur with kernel dimensions

Q 1 = Q 2 = 3.

Influence of each regularization term

We present some numerical and visual results considering different adjustments of ϑ and κ. This experiment allows us to illustrate the influence of the wavelet regularization and the total variation one. In the images displayed in Figure 2, it can be observed the artefacts related to the wavelet regularization, the staircase effects which are typical of the total variation penalization and also the advantage of using a hybrid regularization.

Comparison with existing methods

We now aim at comparing the proposed algorithm with existing methods such as the regularized Expectation-Maximization proposed by Byrne [START_REF] Byrne | Iterative image reconstruction algorithms based on crossentropy minimization[END_REF], the Anscombe approximation method proposed by Dupé et al. [START_REF] Dupé | Image deconvolution under Poisson noise using sparse representations and proximal thresholding iteration[END_REF] and the quadratic extension proposed in [START_REF] Chaux | Nested iterative algorithms for convex constrained image recovery problems[END_REF] (for optimal parameter values). The results are given in Table 2 for "Sebal" image and for different noise intensity factor α. Whatever the scaling parameter chosen, the proposed approach always gives better results in terms of Signal Noise Ratio (SNR). In addition, the improvement in visual quality obtained by adding the total variation penalization to the wavelet regularization is illustrated in Figure 3. α 0.01 0.05 0.1 0.5 EM-Reg. [START_REF] Byrne | Iterative image reconstruction algorithms based on crossentropy minimization[END_REF] 4.12 6.87 8.12 11.1 Anscombe -DTT [START_REF] Dupé | Image deconvolution under Poisson noise using sparse representations and proximal thresholding iteration[END_REF] 7.93 10.8 11.

CONCLUSION

A new approach to restore data degraded by a convolution and Poisson noise has been proposed. The main advan-tages of the method are (i) to deal directly with Kullback-Leibler divergence (without requiring any approximation);

(ii) to permit the use of sophisticated regularization terms, e.g. one promoting sparsity in a wavelet frame and a total variation penalization. Numerical and visual results demonstrate the effectiveness of the proposed approach.

  Therefore, the partition cardinality I introduced in Section 3.2 is such that I 1 = I 2 = 4. C is here defined as (F * ) -1 C * with C * = [0, 255] N . A tight frame version of the dual-tree transform (DTT) proposed in [21] (ν = 2) using Symlets of length 6 has been employed over 3 resolution levels. Strictly convex non-differentiable potential functions are chosen, such that φ k = ω k |.| p k + χ k |.| where (ω k , χ k ) ∈ ]0, +∞[ 2 and p k ∈ {4/3, 3/2, 2}.

Figure 1 :

 1 Figure 1: "Peppers" image (left) and "Sebal" image (right).

Degraded, α = 0. 1 (Figure 2 :

 12 Figure 2: Some results for "Peppers" image.

Figure 3 :

 3 Figure 3: Some visual results for "Sebal" image with α = 0.1.

Table 1

 1 provides quantitative results allowing us to evaluate the impact of the adjustement of the regularization factors ϑ and κ. κ = 0.04 κ = 0.02 κ = 0.01 κ = 0.005

	ϑ = 1 ϑ = 0.2 ϑ = 0.1 ϑ = 0.05 ϑ = 0.01	20.6 21.4 21.6 21.8 21.8	21.6 22.6 22.8 23.0 22.8	21.9 22.9 23.1 23.0 22.1	21.9 22.5 22.1 21.4 19.1

Table 1 :

 1 SNR for "Peppers" image with α = 0.1. Iteration number (IN) lies between 100 and 200.


C Vignat

J.-F Bercher

Geometric aspects of the non-extensive statistical theory

Keywords: nonextensive theory, Tsallis entropy, Rényi entropy PACS: 02.50.Cw

The family of Tsallis entropies was introduced by Tsallis in 1988. The Shannon entropy belongs to this family as the limit case q → 1. The canonical distributions in R n that maximize this entropy under a covariance constraint are easily derived as Student-t (q < 1) and Student-r (q > 1) multivariate distributions. A nice geometrical result about these Student-r distributions is that they are marginal of uniform distributions on a sphere of larger dimension d with the relationship p = n + 2 + 2 q-1 . As q → 1, we recover the famous Poincaré's observation according to which a Gaussian vector can be viewed as the projection of a vector uniformly distributed on the infinite dimensional sphere. A related property in the case q < 1 is also available. Often associated to Rényi-Tsallis entropies is the notion of escort distributions. We provide here a geometric interpretation of these distributions. Another result concerns a universal system in physics, the harmonic oscillator: in the usual quantum context, the waveform of the nth state of the harmonic oscillator is a Gaussian waveform multiplied by the degree n Hermite polynomial. We show, starting from recent results by Carinena et al., that the quantum harmonic oscillator on spaces with constant curvature is described by maximal Tsallis entropy waveforms multiplied by the extended Hermite polynomials derived from this measure. This gives a neat interpretation of the non-extensive parameter q in terms of the curvature of the space the oscillator evolves on; as q → 1, the curvature of the space goes to 0 and we recover the classical harmonic oscillator in R 3 .

THE FAMILY OF EXTENDED ENTROPIES

The family of Tsallis entropies [START_REF] Tsallis | Possible generalization of Boltzmann-Gibbs statistics[END_REF] of a density probability f X is defined by

S q (X ) = 1 1 -q R n f X -f q X , q = 1
where the real positive number q is called the nonextensivity parameter. By l'Hospital's rule we recover the Shannon entropy lim q→1 S q (X ) = -

R n f X log f X
which is thus a limit case. The family of Rényi entropies [START_REF] Rényi | On measures of entropy and information[END_REF] H q (X ) =

1 1 -q log R n f q X
is another set of generalized entropies which contains again the Shannon entropy as the limit case q = 1. The Rényi and Tsallis have the same extremal distributions: since we will focus on these maximal entropy distributions, we will consider indifferently any of these two families. In a classical scheme, we study a complex system whose probability density may evolve according to a complicated law -a differential equation for example. We wish to obtain the equilibrium density, that describes the system after relaxation. A general principle states that for many such systems, the Shannon entropy of the system is increasing

dH 1 (t) dt ≥ 0.
Since moreover the system is often subject to a constraint -for example a fixed value of its energy -its equilibrium probability density can be obtained as the one that maximizes the Shannon entropy under this energy constraint.

For some systems, this H theorem does not hold, but it may be shown that a H q theorem holds [START_REF] Lima | Nonextensive Thermostatistics and the H Theorem[END_REF]: under the same energy constraint, there exists a value of the parameter q such that dH q (t) dt ≥ 0, so that the equilibrium probability density is now the maximum Rényi -or equivalently Tsallis -entropy under the same constraint, hence our interest in the characterization of these distributions.

Let us begin with the explicit expression of the maximum Rényi entropy distributions under a covariance constraint. We remark first that we may assume all variables centered since, with an infinite support, the entropy in invariant by any shift of the variable.

Theorem 1 The maximum Rényi entropy distributions under the covariance constraint E[X X t ] = K

where K is an n × n positive definite matrix are:

• if q > 1, the n-variate Student r-distributions f X;q (X ) = A q 1 -X t Σ -1 X 1 q-1 X ∈ Z ∈ R n |X Σ -1 X ≤ 1 0 else (1) with A q = Γ( p 2 ) Γ( p-n 2 )|πΣ| 1 2
, p = 2q q-1 + n and Σ = pK.

• if n n+2 < q < 1, the n-variate Student t-distributions f X;q (X ) = A q 1 + X t Λ -1 X 1 q-1 ; x ∈ R n (2) with A q = Γ( m+n 2 ) Γ( m 2 )|πΛ| 1 2
; m = 2 1-qn and Λ = (m -2) K.

For obvious reasons, these distributions will be called q-Gaussian distributions, remarking that as q → 1, they converge to the classical Gaussian distribution.

A GEOMETRIC CHARACTERIZATION OF THE MAXIMUM ENTROPY DENSITIES

the case q > 1 A noticeable property of the above distributions is that they can be simply obtained from a universal one, namely the uniform distribution on the sphere, after some simple geometric transformation. We first recall that a random vector U p is uniformly distributed on the sphere S p-1 in R p if it writes

U p = N p N p
where N p is any orthogonally invariant random vector. For example, N p can be chosen Gaussian with identity covariance matrix. We have the following Theorem 2 If U p is uniformly distributed on the sphere S p-1 in R p then any marginal vector X ∈ R n of U p with n < p has a q-Gaussian distribution ( 1) with K = I and q > 1 provided that pn = 2q q -1 .

the case q < 1 An extension of the preceding result can be given in the case q < 1 as follows.

Theorem 3 If the vector X = OM belongs to the unit ball B k , then the point M can be considered as the orthogonal projection on B k of the point P ∈ S k+1 (see Figure 1). The intersection of the line OM with the hyperplane H k+1 = Z ∈ R k+1 |Z k+1 = 1 defines a unique point N such that the vector Y = [ON 1 , . . ., ON k ] t follows a Tsallis distribution with nonextensivity parameter q = n-2 n .

FIGURE 1. the gnomonic projection

We note that the geometric operation that associates the point M to the point N is called a gnomonic projection.

the limit q → 1 case

It is a well-known result that the maximum Shannon entropy n-variate distribution under covariance constraint is the Gaussian distribution

f X;1 (X ) = 1 |2πK| 1 2 exp - X T K -1 X 2 .
This limit is obviously recovered in the expressions (1) and ( 2). However, we may wonder how the geometric approach explicited above extends to this limit case. The answer is given by Poincaré's observation.

Theorem 4 (Poincaré's observation) If U p is a random vector uniformly distributed on the sphere S p-1 and if X n with n < p is any vector of U p then, as p → +∞, the distribution of X n (with fixed n) converges to the Gaussian n-variate distribution with identity covariance matrix.

Our geometric characterization can thus be rephrased as follows: under unit covariance constraint, the maximum Shannon entropy distributions are the marginals of the uniform distribution on an infinite dimensional sphere, while the maximum Rényi entropy distributions are the marginals of the uniform distribution of a finite dimensional sphere.

THE ESCORT DISTRIBUTION

In non-extensive statistics, the escort distribution

f q = f q f q
plays an important role. But its probabilistic meaning is often overlooked; we give here an interpretation of this distribution following a result by Kullback [4,p.40]. Suppose that we have two models for some observed data normalized to a range, say[-1, +1]; the first model is that they follow a distribution f , and the second is that they follow the most natural one, i.e. the uniform distribution U on the interval [-1, +1] . In the absence of other informations, let us decide that the best choice f * for the data is a trade-off between these two distributions. To measure this trade-off, we need a divergence, like the Kullback-Leibler divergence, so that a reasonable way to state our problem is to look for the distribution f * such that

f * = arg min g D (g f ) such that D (g f ) = αD (g U ) (3)
which expresses the fact that f * should be as close as possible to f while remaining close to U in a positive ratio α as well. The solution of the problem (3) is given by the escort distribution f * = f q where the nonextensivity parameter q is related to the ratio α and to the Lagrange parameter λ of the problem by q = 1+λ 1+λ -αλ . Two special cases are (i) α = 0 which obviously implies q = 1 and f * = f , and (ii) α = +∞ for which q = 0 and f * = U. This result shows that the best trade-off (in the sense defined above) between a density f with support [-1, +1] and its uniform counterpart is exactly the escort distribution.

GENERALIZED ENTROPIES AND THE QUANTUM HARMONIC OSCILLATOR

The quantum harmonic oscillator in the plane is the solution of the Schrödinger equation associated with a quadratic potential

Hψ m,n = E m,n ψ m,n
and with Hamiltonian

H = - h2 2m 0 ∆ + 1 2 m 0 ω 2 x 2 + y 2 .
With energy unit hω and length unit h m 0 ω , the solutions are

ψ m,n (x, y) = H m (x) H n (y) e -x 2 +y 2 2 ; E m,n = m + n + 1
and the associated probability density is

f m,n (x, y) = H 2 m (x) e -x 2 H 2 n (y) e -y 2
We remark that the ground state follows the bivariate Gaussian distribution and that the other states are described by the Hermite polynomials, which are orthogonal with respect to the Gaussian measure. These polynomials are described by their Rodrigues formula H n (X ) = (-1) n e X 2 d n dX n e -X 2 .

In a recent paper [START_REF] Cariñena | The quantum harmonic oscillator on the sphere and the hyperbolic plane[END_REF], Cariñena et al. extend these results to the case where the oscillator evolves in a 3-dimensional space of constant curvature κ: the plane when κ = 0, the sphere S 2 when κ > 0 and the hyperbolic plane when κ < 0. Their results are as follows:

Theorem 5 The harmonic oscillator on a space of constant negative curvature

-N (N > 0) is ψ N m,n (z, y) = h N-m-1 2 n (y) h N m (z) ; z = x 1 + y 2 N with h N n (x) = 1 + x 2 N -N 2 -1 4 H N n (x)
and H N n (x) are orthogonal polynomials given by the Rodrigues formula

H N n (x) = (-1) n 1 + x 2 N N+ 1 2 d n dx n 1 + x 2 N n-N-1 2 .
The harmonic oscillator on a space of constant positive curvature ν > 0 is

ψ ν m,n (Z,Y ) = c ν+n+ 1 2 m (Y ) c ν n (Z) ; Z = X 1 -Y 2 ν with c ν n (X ) = 1 -X 2 ν 2 -1 4 C ν n (X ) and C ν n (X )
are the Gegenbauer orthogonal polynomials given by the Rodrigues formula

C ν n (X ) = (-2) n n! Γ (n + ν) Γ (n + 2ν) Γ (ν) Γ (2n + 2ν) 1 -X 2 1 2 -ν d n dX n 1 -X 2 n+ν-1 2 .
From this result, we observe that the ground states in both -negative and positive curvature -cases are exactly q-Gaussian distributions as defined by ( 1) and ( 2) respectively. As a conclusion, the Tsallis entropy extends naturally the role of the Shannon entropy for the harmonic oscillator to the case where the underlying space has constant curvature; moreover, the nonextensivity parameter q can be related explicitely to the curvature in such a way that as q → 1, the curvature equals 0; for more details, see [START_REF] Vignat | Cariñena polynomials are Jacobi polynomials[END_REF].

CONCLUSION

We have shown that a geometric approach can be adopted for some of the aspects of the nonextensive statistical theory: (i) the maximum entropy distributions can be obtained in certain cases as geometric projections of the uniform distribution on the sphere (ii) the escort distribution can be considered as a geometric trade-off in the Kullback Leibler divergence sense and (iii) the extension of the quantum harmonic oscillator to constant curvature spaces involves naturally the nonextensive entropy. More results of this kind are under investigation by the authors in order to enhance the understanding of the nonextensive theory.