Satellite Image Deblurring Using Complex Wavelet Packets.
The deconvolution of blurred and noisy satellite images is an ill-posed inverse problem. Direct inversion leads to unacceptable noise amplification. Usually the problem is regularized during the inversion process. Recently, new approaches have been proposed, in which a rough deconvolution is followed by noise filtering in the wavelet transform domain. Herein, we have developed this second solution, by thresholding the coefficients of a new complex wavelet packet transform&semi; all the parameters are automatically estimated. The use of complex wavelet packets enables translational invariance and improves directional selectivity, while remaining of complexity O(N). A new hybrid thresholding technique leads to high quality results, which exhibit both correctly restored textures and a high SNR in homogeneous areas. Compared to previous algorithms, the proposed method is faster, rotationally invariant and better takes into account the directions of the details and textures of the image, improving restoration. The images deconvolved in this way can be used as they are (the restoration step proposed here can be inserted directly in the acquisition chain), and they can also provide a starting point for an adaptive regularization method, enabling one to obtain sharper edges.
Introduction
The problem presented here is the reconstruction of a satellite image from blurred and
noisy data.
The degradation model is represented by the equation :
Y is the observed data and X the original image. N is additive noise and is assumed
to be Gaussian, white and stationary. The # represents a circular convolution. The Point
Spread Function (PSF) h is positive, and possesses the Shannon property. We deal with
a real satellite image deblurring problem, proposed by the French Space Agency (CNES).
This problem is part of a simulation of the future SPOT 5 satellite. The noise standard
deviation # and the PSF h are assumed known.
The deconvolution problem is ill-posed because of the noise, which contaminates the
data. The inversion process strongly amplifies the noise if no regularization is done.
Thus, we have to deconvolve the observed image while recovering the details, but without
amplifying the noise.
The process used to estimate X from the degraded data Y must preserve textures, to
enable the results to be visually correct. Moreover, the noise should remain small in
homogeneous regions. Many methods have been proposed for regularizing this problem
by introducing a priori constraints on the solution [4], [9], [23]. However, most of them do
not preserve textures, since these textures are not taken into account in the regularizing
model. To achieve a better deconvolution, other authors [1], [19], [21] prefer to use a
wavelet-based regularizing function. But all these approaches have the drawback that
they are iterative, which means that they are time consuming and not always appropriate
for deconvolving satellite data, where the images can be very large.
A few authors, such as Donoho et al. [6], and Mallat and Kalifa [16], proposed denoising
the image after a deconvolution without regularization. The images are represented using
a wavelet or wavelet packet basis, and the denoising process is done in this basis. This
method is not iterative and provides a very fast implementation.
A simple inversion of the observation equation (1) in frequency space gives an unacceptably
noisy solution. To denoise it, a compact representation has to be chosen, in order
to separate the signal from the noise as well as possible. A representation is said to be
compact if it approximates the signal with a small number of parameters, which can be
the coe#cients of the decomposition in a given basis.
The noise amplified by the deconvolution process is colored. Furthermore, the coe#-
cients of this noise are not independent in the wavelet basis. Thus, the basis must adapt
to the covariance properties of the noise. The covariance should be nearly diagonal [15]
w.r.t. the basis, to decorrelate the noise coe#cients as much as possible. The Fourier basis
achieves such a diagonalization, but the energy of the signal is not concentrated over a
small number of coe#cients (the basis vectors are not spatially localized), so the Fourier
transform is not suitable for any thresholding method.
A good compromise is to use a wavelet packet basis [5], since it nearly realizes the
two following essential conditions, i.e. the signal representation is sparse, and the noise
covariance operator is nearly diagonalized [15].
Many types of wavelets transforms can be used to construct a packet basis ; they exhibit
di#erent properties depending on their spatial or frequency localization, and on their
separability w.r.t. rows and columns. Decimated real wavelet transforms are e#cient
for satellite image deconvolution but produce artefacts since the transform is not shift
invariant. To avoid these artefacts, the resulting image has to be averaged over all possible
integer translations. It is also possible to use shift invariant transforms, but the redundancy
is generally very high and depends on the depth of the transform. Anyway, these two
techniques have the major drawback of slowing down the algorithm.
The main motivation of our work is to solve this problem in a computationally e#cient
manner. There is a way to enable translation invariance without much loss of computational
time, by using complex wavelets [17], [18]. Such wavelets also provide a better
restoration by separating 6 directions, while real separable wavelets only take into account
two directions. In order to achieve the necessary near-diagonalization of the deconvolved
noise covariance, we have implemented a complex wavelet packet algorithm.
Our essential contributions are the following
1. We have designed a new transform, the complex wavelet packet transform, which has
better directional selectivity than the complex wavelet transform, while exhibiting the
same shift and rotational invariance properties.
2. The proposed algorithm is fully automatic, since it is based on a Bayesian approach,
where all the necessary parameters are estimated by Maximum Likelihood.
3. We have proposed a new hybrid technique, consisting of combining two di#erent meth-
ods, regularization and wavelet thresholding, to obtain optimal deconvolution results.
4. It performs the inversion much faster than shift invariant real transforms and reconstructs
features of various orientations better.
The paper is organized as follows. First, in section II, we detail how to compute the
proposed complex wavelet packet transform and the properties of this new transform.
Then, in section III, we present the Bayesian thresholding framework used to estimate the
unknown coe#cients from the observed data, and explain how to compute the variance
of the deconvolved noise. Sections III-A to III-C are devoted to the presentation of the
di#erent prior models put on the wavelet coe#cients - homogeneous, noninformative and
inhomogeneous priors. In section III-C.2, we detail the hybrid technique used to estimate
the adaptive parameters of the latter. This model is used in the algorithm proposed in
section IV. Finally, we present, in section V, a comparison with classical algorithms used
for satellite image deconvolution, to demonstrate the superiority of the proposed method.
II. Complex wavelet packets
To build a complex wavelet transform, Kingsbury [17] has developed a quad-tree al-
gorithm, by noting that an approximate shift invariance can be obtained with a real
biorthogonal transform by doubling the sampling rate at each scale. This is achieved by
computing 4 parallel wavelet trees, which are di#erently subsampled. Thus, the redundancy
is limited to 4, compared to real shift invariant transforms.
The shift invariance is perfect at level 1, and approximately achieved beyond this level :
the transform algorithm is designed to optimize the translation invariance. Therefore, it
involves two pairs of biorthogonal filters, odd, h o and g o , and even, h e and g e . At level
simply a non-decimated wavelet transform (using h
are re-ordered into 4 interleaved images by using their parity. This defines the 4 trees
B, C and D. For j > 1, each tree is processed separately, as a real transform, with
a combination of odd and even filters depending on each tree. The transform is achieved
by a fast filter bank technique, of complexity O(N ).
We have extended the original transform by applying the filters h and g on the detail
subbands, thus defining a complex wavelet packet (CWP) transform [13], [14]. This
new transform exhibits the same invariance properties as the original complex wavelet
transform. The tree corresponding to this transform is given in Fig. 1.
The filter bank used for the decomposition is illustrated by Fig. 2, on which the subbands
are indexed by (p, q) for each tree T . The impulse responses, shown in Fig. 3, and
the related partitioning of the frequency space given in Fig. 4, demonstrate the ability
to separate up to 26 directions for the chosen decomposition tree. Compared to real
separable transforms, which only define two directions (rows and columns), it provides
near rotational invariance and gives a selectivity which better represents strongly oriented
textures (thus separating them better from the noise). Furthermore, compared to the
original complex wavelet transform, which only separates 6 directions, the directional
selectivity is higher. Instead of dividing an approximation space which does not define any
new orientation, the wavelet packet decomposition processes the detail subbands, which
are strongly oriented. Each detail subband at level 1 isolates an area of the frequency
space defined by a mean direction and a dispersion, enabling one to select a range of
directions around an orientation #. If the subband is decomposed into 4 new subbands,
it means that the corresponding frequency area splits into 4 new areas, which can define
new orientations, as shown in Fig. 4.
This is not really a complex transform, since it is not based on a continuous complex
mother wavelet. Nevertheless, the quad-tree transform has in practice the same properties
as a complex transform w.r.t. shifting of the input image, and it is perfectly invertable
and computationally e#cient. The complex coe#cients are obtained by combining the
di#erent trees together. If we index the subbands by k, the detail subbands d j,k of the
parallel trees A, B, C and D are combined to form complex subbands z j,k
and z j,k
- , by a
linear transform (represented by a matrix M ), in the following way :
z j,k
A - d j,k
z j,k
Thresholding the magnitudes |z - | without modifying the phase enables us to define a
nearly shift invariant filtering method.
The reconstruction is done in each tree independently, by using the dual filters, and the
results of the 4 trees are averaged to obtain a 0 to ensure the symmetry between the trees,
thus enabling the desired shift invariance.
III. Bayesian thresholding
Let us denote X the deconvolved image without regularization. For each subband k,
the variables x and # denote one of the CWP coe#cients corresponding respectively to the
observation Y and the original image X . We suppose that H is invertible (i.e. it has no
zeros in the Fourier space). Since the noise N is white and Gaussian, Equ. (1) multiplied
by H -1 gives, in the CWP domain :
To obtain the expression of the noise coe#cients n, let us recall Equ. (2). Each complex
coe#cient is obtained by summing or subtracting the coe#cients of the 4 trees A, B, C,
D. If we compute the covariance between the real and imaginary parts, we find (for the
coe#cients z are the
noise coe#cients for each tree. By symmetry assumptions, this covariance is null, since all
covariances E[n T n T # ] between di#erent trees T and T # are equal. Then, the distribution
of the noise is defined as a joint distribution of (n r ,
We assume that the coe#cients # are independent in a given subband, between di#erent
subbands of a given scale, and also between scales. This is an approximation which
enables a fast thresholding technique: we will not handle here possible correlations between
subbands or neighbour coe#cients. The covariance matrix of the noise is supposed to be
nearly diagonal in the chosen basis (see Fig. 6 and Fig. 7), so we consider that the noise
variables are also independent in the wavelet packet basis.
We assume that the noise variance is constant in each subband k. We compute # k by
considering the undecimated transform of the noise N , which is performed by a linear
operator (convolution with the impulse response w k , obtained by an inverse CWPT
of a Dirac). F denotes the Fourier transform.
We estimate the unknown coe#cients # within a Bayesian framework [22]. We have
demonstrated in [14] that this approach provides slightly better results than the Minimax
risk calculus on real satellite data. To compute the MAP estimate of #, we use Bayes
law to calculate the expression of the posterior probability :
where P (x | #) is given by the distribution of the noise :
A. Homogeneous prior model
Generalized Gaussian distributions have been used to model real wavelet coe#cients
[20], [22]. We also propose to use them to model the CWP subbands. We have the prior
probability of # :
where # k is a prior parameter and p k is an exponent. However, we assume that the variables
# are independent within a given subband, and also between the di#erent subbands. As
shown by Fig. 5, the complex density is a bidimensional function which only depends on
the magnitude (it exhibits a radial symmetry). It is generally not separable for p #= 2. We
set a given density function on the magnitude, while the phase is uniformly distributed in
The exponent p k can be set to a fixed value, the same for all subbands, to simplify the
computation. Indeed, if this parameter is not specified it must be estimated. This is quite
complex and is not justified by the improvement of the results.
On large size images, the parameter # k can be computed e#ciently by various methods,
such as Maximum Likelihood for example. We choose to estimate it automatically from
the histogram of a given subband.
Then we obtain the following expression of the MAP, using the prior law
It is possible to demonstrate that it is equivalent to apply a thresholding operator # p k
to the magnitude of each coe#cient [14] :
For complex wavelets, this leads to the following equations :
and
This shows that # is obtained from x by keeping the phase, and thresholding the magnitude.
It reduces to simple soft thresholding if
elsewhere. Thus, we have :
This classical filtering method naturally arises from the Bayesian approach. The resulting
thresholding functions are smoother for is the exponent of the prior law (7)),
and they become linear if 2. If p k < 1, which is more realistic for satellite images, the
functions become discontinuous. Then, the thresholding functions # are numerically
computed by solving equation (10) w.r.t. Fig. 8 shows the behaviour of these functions
for di#erent values of p k . Experimental studies have shown that provides an
e#cient model for satellite images [14].
B. Noninformative Je#rey's prior
It is possible to use a di#erent approach for subband modeling. This approach has
been proposed in [8] within a wavelet based image denoising framework. It is based on
the following assumption : the inference procedure should be invariant under changes of
amplitude and scale. It means that the prior probability law of # must keep the same
behaviour even if it is rescaled, by setting for example this is not the case
for classical Gaussian or Laplacian models, the author uses the following prior, which is
called a noninformative prior
This corresponds to an extremely heavy-tailed distribution, which approximately describes
the original wavelet coe#cients #. Unfortunately it is improper : the resulting posterior
density function is not integrable.
Therefore, an alternative to the fully Bayesian framework is chosen. It consists of
treating the real or imaginary part of each coe#cient as a zero-mean Gaussian variable, of
variance s 2 . This defines an adaptive model, since each coe#cient has a di#erent variance.
Then, this variance is supposed to follow Je#rey's hyperprior distribution, i.e. we have :
This is equivalent to P (# 1/|#| within a Bayesian context. Thus, the model remains
homogeneous, even though an adaptive Gaussian model is used intermediately to address
the problem of the improper distribution.
The estimation is then performed in two steps :
. estimate the variance - s 2 by using the
| x)
. estimate the unknown coe#cient by using the
To express P (s 2
| x), we need P both signal and noise are
Gaussian, of respective variances s 2 and # 2
k , we have P
We also have the hyperprior
(13). Then we
|
which is equivalent to s
The MAP estimate -
for the inhomogeneous Gaussian model gives the
following expression (see the next section for a complete proof) :
By combining equations (14) and (15) we obtain the following estimate, which we call the
noninformative thresholding function # J :
The advantage of such a method is that there is no need for parameter estimation,
since there is no parameter. However, there is a drawback, which is especially visible in
homogeneous areas : the residual noise is quite apparent, its variance being higher than
with the previously presented homogeneous model. This probably comes from the lack of
robustness of the estimation method. We can remark that if we remove the prior law of
the estimation of s 2 is done by the MLE, which gives a function like
but with a threshold 2# 2 instead of 4# 2 . This is insu#cient since the magnitude of
the noise has a variance equal to 2# 2 . Thus, the hyperprior makes the estimation more
robust, by doubling the threshold value. It is still not su#cient to remove noise peaks
in large constant areas. Hereafter, we detail a more robust method to filter the complex
wavelet packet subbands.
C. Inhomogeneous prior model
C.1 The model
The Generalized Gaussian model presented previously has been chosen because it seems
to match correctly the original coe#cient distribution, which is heavy-tailed. Another
possible way to capture this property is to define an inhomogeneous model, which adapts
to the local characteristics of the subbands. Some approaches to spatial adaptivity in the
real wavelet domain can be found in the literature, for example [3].
To simplify, let us choose a Gaussian model. The variance parameter is di#erent for
each coe#cient, which enables us to di#erentiate edges or textures, which have a high
intensity, from the homogeneous areas which generally correspond to very low values of
the coe#cients.
Since the parameters can be very di#erent from one variable to another, the histogram
of a subband can have a heavy-tailed behaviour, even if the distribution of each variable
is Gaussian.
We denote by s 2 the variance of the real or imaginary part of an original coe#cient #,
as before. We have the prior law :
If the parameters of the prior distribution are known, the unknown variables are estimated
by computing the MAP. This is a fully Bayesian technique (we use the property (5)). Recall
the expression of the noise distribution (6), and combine it with Equ. (17) to obtain :
Therefore, the MAP is given by :
To compute the minimum, di#erentiate w.r.t. the real and the imaginary parts of #. This
gives two equations which are recombined to form an equation with complex numbers :
#r +xr
This gives the inhomogeneous MAP estimate :
C.2 Adaptive parameter estimation
The most di#cult problem in this approach is to estimate the adaptive parameters of
the model. As is shown in [12], the MLE is not robust when applied to incomplete data
(i.e. when the estimation is made on noisy data). Indeed, there are as many parameters as
observed data. But the robustness becomes su#cient when the estimation is made from
the original image X . A good approximation of this image is still su#cient to provide
useful parameter estimates. Here, we use this complete data approach to estimate the
variances s 2 .
Consider Equ. (17). The complete data MLE is defined by - s
Assumimg independence, we obtain :
where the factor 2 comes from the dimensionality of the distribution.
We obviously do not have access to the original coe#cients, which is why we take the
transform coe#cients of an approximate original image instead. Experiments have shown
that a satisfactory approximation is provided by a nonlinear regularizing algorithm, such
as RHEA, detailed in [11]. It essentially consists of a variational method based on #-
functions [4] (minimization of a criterion which penalizes noisy solutions, but preserves
edges) preceded by an automatic parameter estimation step to compute the hyperparameters
of the regularizing model.
This method is certainly not perfect and some residual noise remains. It is visible
in constant areas. However, we filter this noise as well as the deconvolved noise by a
thresholding technique. We choose to use the noninformative threshold of the previous
section, because it does not require any additional parameter estimation.
The proposed algorithm consists of obtaining the desired approximate original image
using RHEA [11], filtering the CWPT of the result using Equ. (16), estimating the adaptive
parameters using the complete data MLE with Equ. (20), and then estimating the
unknown coe#cients by computing the MAP by Equ. (19).
In addition to the computation of the deconvolved noise variance, we also need to
compute the variance 2-# 2
k of the residual noise of the approximate image.
Let us denote by -
# the thresholded transform coe#cients of the approximate original
image. -
# is supposed to be su#ciently exempt from noise and to contain su#cient information
to enable texture and edge recovery. Homogeneous areas and edges are fine, since
we have used an edge-preserving method followed by an e#cient noise thresholding. But
we still have to explain why the coe#cients related to textured areas are su#ciently high
to avoid too strong an attenuation by using Equ. (19) in these areas. The variational
method used does not completely remove the textures, and even if they are visually not
very sharp, they are su#ciently present in the approximate image to enable a correct
reconstruction using the method detailed here.
Finally, if -
# is known, by using Equ. (19) and (20), the estimate for the coe#cient is :
| -
If we compute the expression of the thresholding function which minimizes the Bayesian
risk of the estimator -
#, defined by
denotes an expectation w.r.t. the distribution of x, we find :
It has exactly the same form as Equ. (21), if we take the approximate original coe#cient
# instead of #. Both Bayesian and minimum risk methods lead to the same expression,
which means that the chosen model is good, since the corresponding estimator provides
the minimum risk.
As in the case of Wiener filtering [10], computing the MAP under Gaussian assumptions
(both signal and noise have Gaussian distributions) is equivalent to minimizing the risk
of a linear estimator w.r.t. the attenuation factor. Therefore, the two approaches are
equivalent in the Gaussian case.
IV. The deconvolution algorithm
The adaptive model described by Equ. (17) provides much better results (visually and
w.r.t. SNR) than the homogeneous model described by Equ. (7). Details are better
preserved and constant areas are cleaner. That is why we keep this model for the final
version of the proposed algorithm. We have also compared this scheme with the classical
approach [6] and with minimum risk computation for various models [14]. It
consistently exhibits better results.
The initial deconvolution is made in the Cosine Transform space instead of the Fourier
space to avoid artefacts near the borders of the image. The proposed algorithm is called
COWPATH, for COmplex Wavelet Packet Automatic THresholding, and consists of the
following steps (see Fig.
. DCT (Discrete Cosine Transform) of the observation Y
. Deconvolution : divide by F [h] (in practice, divide by F [h] is small, since
some of the coe#cients of F [h] can be null)
. Inverse DCT of the result, which gives X
. CWP transform of X
. Computation of the approximate original image -
apply the RHEA algorithm [11] on
Y (nonlinear regularization, with automatic parameter estimation)
. CWP transform of -
. Computation of # k using the known h and # (see Equ. (4))
. Computation of - # k (residual noise on the approximate original image) using the known
h and # (see [14] for details)
. Thresholding of the approximate image coe#cients using Je#rey's noninformative prior
and -
. Estimation of the parameters - s of the inhomogeneous Gaussian model (see Equ. (20))
. Coe#cient thresholding by computing the MAP (see Equ. (21))
. Inverse CWP transform, which gives the estimate -
X.
The variance of the residual noise in homogeneous areas, which is needed to denoise the
approximate original image before using this image for parameter estimation, is computed
in the same manner as the variance of the deconvolved noise (i.e. using an equation
similar to Equ. (4)). For this computation, we assume that constant regions correspond
to a quadratic regularization. Then it is possible to use sums in the Fourier space. We
refer to [14] for more details.
The main characteristic of this algorithm is the use of two di#erent methods, regularization
and wavelet thresholding, to obtain a hybrid technique whose results are better than
the results of a single regularization method or a wavelet thresholding. The more decorrelated
the deconvolved noise and the residual noise of the approximate original image, the
higher the quality of the deconvolved image.
It is possible to replace the nonquadratic regularizing model of the RHEA algorithm by
a simple quadratic model. The advantages are to enable a single step deconvolution in the
frequency space and to make the parameter estimation step deterministic and fast (in the
nonquadratic case we need a MCMC method [11]). The edges are not as sharp as with
the nonquadratic model, but this approximate image is su#ciently accurate to provide
a correct estimation of the inhomogeneous parameters of the subband model (the SNR
di#erence between the original and accelerated algorithms is about 0.1 dB for the SPOT
5 simulation image shown in Fig. 13).
The complexity of the algorithm is 115 log (operations per pixel) (1980
op/pix for a 512 - 512 image). If we use the accelerated version, based on the previously
described quadratic model, it is 15 log
image, which represents about 4 s on a Pentium II 400 MHz machine).
V. Satellite image deblurring
A. simulation
Fig. 10 shows a 128x128 area extracted from the original image of N- mes (SPOT 5
simulation at 2.5 m resolution, provided by the French Space Agency (CNES)). Fig. 11
shows the PSF and Fig. 12 shows the observed image Y
Fig. 13 shows the image deconvolved with the proposed algorithm.
B. Comparison with di#erent methods
1. Quadratic regularization [23]. This is nearly equivalent to the parametric Wiener filter
which gives the same results. It is also equivalent to isotropic di#usion. The edges
are filtered as well as the noise, as seen on Fig. 14. It is therefore impossible to obtain
sharp details and noisefree homogeneous areas at the same time. Thus, the SNR remains
low (about 19.7 dB) because of insu#cient noise removal in these areas.
2. Nonquadratic regularization [4], [9]. The resulting image of the RHEA algorithm [11]
exhibits sharp edges, compared to the result shown previously. However, some noise
remains in homogeneous regions and textures are attenuated. The SNR is 22.0 dB. This
result is used as the approximate original image and is illustrated by Fig. 15.
3. Real wavelet packets [15]. The proposed complex wavelet packet transform is more than
two times faster than the shift invariant real wavelet packet transform (based on Symmlet-4
wavelets) and is much more directionally selective. Real wavelet packet thresholding gives
a SNR equal to 21.8 dB. See Fig. 16 for an illustration.
4. The proposed method. This is faster than the other methods, and provides the highest
SNR (22.2 dB). The textures and the oriented features are sharp and regular, while the
homogeneous regions remain noise free, as seen on Fig. 13.
VI. Conclusion
We have proposed a new complex wavelet packet transform to make an e#cient satellite
image deconvolution algorithm. This transform exhibits better directional and shift
invariance properties than real wavelet packet transforms, for a lower computational cost.
The proposed deconvolution method is superior to other competing algorithms on satellite
images : it is faster, more accurate and fully automatic. The essential novelty of
the proposed algorithm consists of an hybrid approach, in which two radically di#erent
methods are combined to produce a deconvolved image of higher quality than the result
of each method individually. Finally, if a quadratic model is used, the speed is greatly
increased, which opens the path to real time processing of image sequences or to onboard
image processing in satellites by using specialized chips.
Furthermore, this new type of approach can be extended and the results can be improved
to handle more di#cult cases (di#erent types of blur and higher noise variance). It is
possible to take into account various levels of noise and di#erent convolution kernels.
Indeed, the wavelet packet decomposition tree is not unique and should be adapted to the
noise statistics.
The case of noninvertible blur, such as motion blur, forbids the use of nonregularized
inversion in the Fourier domain. Therefore, in this case, a regularized deconvolution should
replace the rough inversion used in the proposed method.
An improvement of the adaptive Gaussian model could also be provided by choosing a
more accurate distribution, to capture the heavy-tailed distribution of the wavelet packet
coe#cients. For example, adaptive Generalized Gaussian models could be investigated.
Furthermore, including dependence between di#erent scales by means of multiscale hidden
Markov trees could perhaps enable better separation of the small features from the deconvolved
noise. These types of models would certainly better model edges which propagate
across scales and enable their reconstruction more e#ciently.
VII.

Acknowledgements

The authors would like to thank J-er-ome Kalifa (from CMAPX, at Ecole Polytechnique)
for interesting discussions and Nick Kingsbury (from the Signal Processing Group, Dept.
of Eng., University of Cambridge) for the complex wavelets source code and collaboration,
Peter de Rivaz (same institution) for his kind remarks, and the French Space Agency
(CNES) for providing the image of N- mes (SPOT 5 simulation).


--R

Wavelet domain image restoration using edge preserving prior models.
Bayesian Theory.
Spatial adaptive wavelet thresholding for image denoising.

Wavelet analysis and signal processing
Denoising by soft thresholding.
spatial adaptation via wavelet shrinkage.
Bayesian wavelet-based image estimation using noninformative priors
Stochastic Relaxation
Image restoration by the method of least squares.




Restauration minimax et d-econvolution dans une base d'ondelettes miroirs
Wavelet packet deconvolutions.
The dual-tree complex wavelet transform: a new e#cient tool for image restoration and enhancement
"Wavelets: the key to intermittent information?"

A theory for multiresolution signal decomposition: the wavelet representation.
A wavelet regularization method for di
Bayesian inference in wavelet based methods
Regularization of incorrectly posed problems.
--TR
complex wavelet packets;satellite and aerial images;deblurring;bayesian estimation