Regularized Richardson-Lucy Algorithm for Sparse Reconstruction of Poissonian Images
Restoration of digital images from their degraded measurements has always been a problem of great theoretical and practical importance in numerous applications of imaging sciences. A specific solution to the problem of image restoration is generally determined by the nature of degradation phenomenon as well as by the statistical properties of measurement noises. The present study is concerned with the case in which the images of interest are corrupted by convolutional blurs and Poisson noises. To deal with such problems, there exists a range of solution methods which are based on the principles originating from the fixed-point algorithm of Richardson and Lucy (RL). In this paper, we provide conceptual and experimental proof that such methods tend to converge to sparse solutions, which makes them applicable only to those images which can be represented by a relatively small number of non-zero samples in the spatial domain. Unfortunately, the set of such images is relatively small, which restricts the applicability of RL-type methods. On the other hand, virtually all practical images admit sparse representations in the domain of a properly designed linear transform. To take advantage of this fact, it is therefore tempting to modify the RL algorithm so as to make it recover representation coefficients, rather than the values of their associated image. Such modification is introduced in this paper. Apart from the generality of its assumptions, the proposed method is also superior to many established reconstruction approaches in terms of estimation accuracy and computational complexity. This and other conclusions of this study are validated through a series of numerical experiments.
The notion of event counts is fundamental in many imaging modalities. Thus, for instance, this notion is used to quantify the reception of gamma photons in nuclear imaging [1, 2, 3, 4] as well as to describe the formation of optical images in charge coupled devices (CCD) [5, 6]. Confocal microscopy , astronomical  and turbulent imaging  are additional examples of applications in which the notion of even counts routinely arises. In all these cases, working with event counts entails using a specific statistical assumption on the nature of acquired images. In particular, the images are assumed to be formed as blurred versions of their original counterparts contaminated by Poisson noise. Consequently, to recover the original images, the combined effect of the blur and noise needs to be reversed. A novel approach to the solution of this classical inverse problem is proposed in this work.
Since in practice, digital images are discrete and finite-dimensional objects, it seems reasonable to represent the original image along with its measured version by real-valued matrices. Moreover, since in Poissonian imaging the values of both and have the interpretation of event counts, the images can be further constrained to be members of the subspace of nonnegative-valued matrices. Formally,
Additionally, let be a linear operator describing the convolution of with a non-negative mask of some size. Then, the formation of can be formally expressed as
where describes the effect of contamination of by Poisson noise. Therefore, as was mentioned before, to recover the original image from , the combined effect of in (2) needs to be inverted.
The current arsenal of image restoration approaches to the solution of (2) is relatively broad. In general, all these methods can be divided into two main groups, where the methods of the first group assume to be identity, whereas the methods of the second group allow it to be a general low-pass filter. Within the first group, a range of available reconstruction methods take advantage of a variance stabilized transformation (VST) , which allows the Poisson noise in (2) to be transformed into approximately white Gaussian noise. Some recent developments in this direction include the works in [11, 12, 13], in which wavelet de-noising has been exploited to reject the “gaussianized” Poisson noise. Conceptually similar ideas have been also advocated in [14, 15] based on more advanced tools of statistical analysis. The framework of Bayesian estimation was exploited in [16, 17, 18, 19, 20], where hierarchical (Markovian) models have been used to describe a priori probabilities of the original images. Although all the above methods can be used to recover an approximation of in (2), their applicability is restricted to the case of weak blurs, where can be closely approximated by the identity operator . Moreover, all these methods depend on VST, whose performance is known to deteriorate considerably in low-count settings .
One of the nowadays classical methods of solving (2) for the case of a general was proposed in the beginning of the 1970s in the seminal works of W. Richardson  and L. Luci . This method – known as the Richardson-Lucy (RL) algorithm – can be classified as a maximum likelihood (ML) estimator. As a general rule, however, ML estimation may result in less accurate and/or stable reconstructions as compared with maximum a posteriori (MAP) estimation methods based on the Bayesian paradigm. Consequently, the RL algorithm has been recently extended under the MAP estimation framework, resulting in a number of regularized solutions to (2) which differ in the way the original image is modelled. Thus, for example, Gaussian models have been exploited in [24, 25], while the algorithm in  is based on the total variation model of . Unfortunately, there are conditions on which the above methods can result in erroneous reconstructions (as will be demonstrated later in the paper).
Despite the relative simplicity of the RL algorithm, it still remains one of the most widespread methods used in the current practice of Poissonian imaging. The main advantage of the RL algorithm stems from its nonlinear nature, allowing one to recover the high-frequency components of the original images which are the most affected by blur. Moreover, a closer look at the analytical properties of the RL algorithm reveals the fact that it has a “built-in” ability to converge to sparse solutions, and, as a result, the application of the algorithm for the restoration of sparse images should be expected to produce particularly useful reconstructions. Unfortunately, most of the real-life images are not sparse in the spatial domain, in which case the RL algorithm may fail to provide useful results. On the other hand, most of such images are sparse in the domain of a properly designed linear transform [28, 29, 30, 31]. Consequently, to extend the applicability of the RL algorithm to a wider range of imagery data, it is tempting to find a way to apply the algorithm in the transformed domain, as opposed to the spatial domain. Accordingly, the main contribution of this paper consists in the introduction of a novel approach to the solution of (2) which exploits the above idea. Moreover, it will be shown via an experimental study that the proposed method outperforms a number of alternative algorithms in terms of normalized mean-square error (NMSE), SSIM quality index , as well as its stability and computational efficiency.
The rest of the paper is organized as follows. Section II provides some necessary details on the RL algorithm. The statistical models and assumptions underpinning the proposed method are detailed in Section III, while Section IV summarizes the main structure of the proposed solution method. The results of both computed-simulated and real-life experiments are presented in Section V. Section VI finalizes the paper with a discussion and conclusions.
2 Richardson-Lucy Approach
To establish a necessary foundation for the development and presentation of the proposed method, a brief overview of the RL algorithm is provided first. Under the image formation model of (2), the likelihood function of the measured image can be defined as
where the samples of are assumed to be mutually independent. Consequently, defining , the ML estimate of is given by
with denoting the standard inner product in . Differentiating with respect to and equating the derivative to zero results in
where denotes an matrix of ones, the fractional line stands for a point-wise division, and denotes the adjoint of . Note that if represents the operation of 2-D convolution with a positive-valued kernel , then represents the convolution with a spatially reversed version of the same kernel . Moreover, if is normalized to satisfy , then both and are obviously equal to . This fact is exploited by the RL algorithm, which recovers an approximation of the original image as a stationary point of the sequence of solutions produced by
where the dot stands for a point-wise multiplication, and is the iteration index.
To better understand the nature of the solution produced by the minimization of (5), consider the following steps. Using the definition of an adjoint operator, the cost functional in (5) can be rewritten as
Moreover, since is assumed to be positive valued, the inner product coincides with the -norm of defined as . As a result, the cost function in (4) can be further redefined as
A closer look into the two terms in (9) leads to a number of remarkable observations. First and foremost, the presence of -norm in the cost function implies that the solution of (4) will tend to be sparse [28, 29, 33, 34, 35]. Therefore, the RL algorithm has a “built-in” ability to converge to sparse solutions. Moreover, such solutions are guaranteed to be positive due to the second term in (9) which works akin to a logarithmic barrier [36, Ch.11]. This is what makes the RL algorithm to favour the reconstructions which are sparse in the spatial domain.
It goes without saying that the assumption of spatially sparse images may be unacceptable for a variety of natural scenes. For this reason, the RL algorithm is known to perform reliably in the case of, e.g., astronomical images , while its application to piecewise smooth images can produce rather disappointing results . On the other hand, sparsity is known to be an extremely useful and liable assumption to use for describing the behaviour of the representation coefficients of natural images in the domain of certain linear transforms. Therefore, to extend the applicability of the RL algorithm to the above case, it is imperative to find a way to apply the algorithm directly to the (sparse) representation coefficients, rather than to their corresponding images. A possible variant of such an approach is detailed in the section that follows.
Before turning to the description of our method, a number of important properties of the RL recursion in (7) are worth to be mentioned. In particular, it can be seen from (7) that the algorithm preserves the positivity of , viz. is guaranteed to be in provided and . Somewhat less trivial, it can also be shown that the algorithm preserves the mean values of the intermediate solutions, i.e. at each iteration . All these properties will play an important role in the proposed method as detailed below.
3 MAP Formulation
The approach proposed in this paper is based on the assumption that the original image in (2) can be sparsely represented in the domain of a linear transformation. In particular, given an overcomplete and dense set of vectors in , we define the synthesis operator to be given by
where denotes the set of indices of the representation coefficients , with . Using this definition of , one can rewrite the model equation (2) as
with being a composition of and , i.e. , and being the representation coefficients of in the domain of . In this new format, the original problem of reconstruction of is replaced by the problem of estimating its representation coefficients .
Needless to say, due to the over-completeness of , the definition of is not unique. This difficultly, however, can be easily overcome by requiring the representation coefficients to be as sparse as possible. A possible way to quantify the sparseness of is by measuring its norm . It is worthwhile noting that this measure of sparseness has been successfully used in numerous fields of science, which include signal modeling , compressed sensing [34, 35], independent component analysis and blind source separation [37, 38], inverse problems [39, 40], signal and image de-noising [28, 30], morphological component analysis  together with its earlier version introduced in .
Using the assumptions above, the problem of recovering the original image can be cast into the framework of MAP estimation, in which case the optimal is found as a solution of the following maximization problem
where and denote the likelihood function and the prior probability of , respectively.
Under the assumption of statistical independence, the likelihood has the form of
One the other hand, congruent to the assumption on minimality of , the representation coefficients are assumed to be i.i.d. random variables that obey a zero-mean Laplacian distribution, viz
where is a scale parameter of the distribution.
with being a regularization parameter.
Care should be exercised when specifying the domain of definition of in (16). Indeed, the likelihood model of (13) interprets the value as the mean of the corresponding random observation . Since, in the case of a Poisson distribution, its first and second moments are equal, the values should be assumed to be non-negative (in accordance with our earlier assumption in (1)). Thus, the domain of the objective is formally given by
which is a convex set. Moreover, as long as and the convolution operator is non-degenerate (albeit, possibly, ill-conditioned), the objective functional is guaranteed to be strictly convex. In this case, admits a unique minimizer in , which can be found by any optimization algorithm which is guaranteed to converge to a stationary point of . Unfortunately, when either is not strictly positive or possesses a non-trivial null-space, the convexity of is not strict, and, as a result, the existence and uniqueness of its global minimizer cannot be a priori guaranteed.
4 Solution via Fixed-Point Iterations
Numerous algorithms approaches can be used to minimize the cost functional in (16). In this paper, we propose a different method based on fixed-point iterations. To this end, we first notice that the first-order optimality condition for (16) has the following form
where denotes the adjoint of and is an extremum of . Note that, in (18), the non-differentiability of the absolute value at zero is resolved by letting . As will be shown below, the latter assumption is rather technical, as it has no impact on the proposed solution.
Finally, by rearranging the terms in (18), one obtains
which, in turn, suggests the following fixed-point iteration algorithm
The iterative procedure of (20) is multiplicative in nature – the fact that has a number of important implications. First, the zero entries of are preserved by the iteration procedure, which means that if for some , then for all . In this respect, letting seems to be a reasonable simplification, since the value of the gradient of at zero does not seem to have any influence on the result of the iterative procedure (20). Second, if the values of were allowed to be of an arbitrary sign, it would be generally impossible to guarantee a monotone convergence of . To overcome this deficiency, in what follows the representation coefficients will be assumed to be non-negative. Moreover, if the representation vectors are constrained to be positive-valued as well, it is straightforward to verify that the iterations in (20) will preserve the positiveness of . It is important to note that the above constraints on the sign of and should not be seen as a limitation, since there exists a body of works (see, for example, ), which describe the construction of positive-valued representation sets that allows the natural scenes to be represented in terms of both sparse and positive coefficients.
Provided the elements of are positive values, the vector will be positive valued as well, i.e. . In this case, (20) can be rewritten in a more compact form as
The iterations in (21) can be initialized with a constant coefficient vector (e.g. ) and executed until the relative change drops below a predefined threshold . The above algorithm will be referred below to as the sparse RL method, or simply SRL.
Comparing the recursions in (21) and (7), one cannot avoid noting how similar the SRL and RL algorithms are. Indeed, replacing and in (21) by and , respectively, yields an update equation identical to that in (7) (up to the element-wise division by in (21)). Moreover, the assumption on the non-negativeness of and made by SRL is parallel to the non-negativeness of and exploited by the RL method. Further similarity between RL and SRL reveals itself in the objective functions of the two reconstruction methods. Specifically, the objective function in (16) can be rewritten as
where and is the -weighted -norm of the non-negative representation coefficients . Thus, one can see that under the substitution of by , by and , the objective functionals (4) and (9) are identical. Yet, while minimizing (9) forces the reconstruction of to have a sparse appearance, minimizing (4) has the same effect on the representation coefficients. As opposed to the former, the latter case allows one to recover much broader classes of practical images, as it is demonstrated by the experimental results below.
5 Experimental Results
5.1 Reference methods and comparison measures
Both one- and two-dimensional data sets have been used to test the performance of the proposed SRL algorithm. In the case of 1-D data, the SRL algorithm was compared against the RL method, while in the case of 2-D data, the comparison was made against the method described in . Note that, just like SRL, the reference method of  belongs to the class of Bayesian estimators. Specifically, the method recovers the original image as a maximizer of its posterior probability, computed under the model of (2) and the total-variation (TV) prior. The maximization can be performed iteratively according to the following update equation
where is a regularization parameter, which has been set to be equal to 0.002 as prescribed in . For the convenience of referencing, the above algorithm will be referred below to as the RLTV method.
To quantitatively compare between the performance of various reconstruction algorithms, two comparison measures have been used. The first, normalized mean-squared error (NMSE), is defined as follows. Let be an original image and be an estimate of . Then, the NMSE can be defined as
with being the Frobenius matrix norm, and being the operator of expectation. In the current study, the expectation was approximated by its corresponding sample mean based on the results of 200 independent trials.
It has been recently argued that the NMSE may not be an optimal comparison measure as long as human visual perception is concerned. For this reason, the structural similarity (SSIM) index proposed by  was employed as well as an alternative performance metric.
5.2 One-dimensional Reconstruction
The first part of our experimental study is concerned with the recovery of piece-wise constant signals of length . For this case, a natural choice would be to define the basis functions to be scaled and spatially shifted versions of a rectangular (Haar) box. Specifically, let be a (dyadic) resolution index, and be defined as
where the normalization factor is used to guarantee for all . Also, let
be the operator of (causal) circular shift by points. Obviously, for each there are a total of non-repetitive shifts. Consistent with the notations standard in wavelet theory, let denote the function (circularly) shifted by points to the right (in which case can be viewed as a rectangular box function supported on ). Moreover, let be a matrix of size , whose columns are formed by functions with and . In this case, given an arbitrary vector of representation coefficients of size , its corresponding signal can be synthesized as .
The dictionary constructed above is severely overcomplete, which has the disadvantage of high computational complexity associated with the computations of and . To improve the computational efficiency, the index was restricted to four resolution levels, viz. , and . As a result, the size of was equal to for . The coefficients used for the synthesis of test signals had about % of non-zero entries, drawn from a uniform distribution. The blurring artifact was simulated by convolving the test images with a Gaussian kernel whose -3 dB cut-off frequency was set to be equal to . The maxima of the resulting signals were set to two different values, namely 256 and 32, in order to simulate high- and low-count detection scenarios, respectively. As a final step, the blurred signals were contaminated by Poisson noise (see the second subplots in Fig. 1 and Fig. 2 for typical examples of data signals).
Next, the RL and SRL algorithms were applied to the synthesized data signals with . Unfortunately, no monotonous convergence in NMSE could be achieved in the case of the RL algorithm. The steady-sate estimates obtained using this classical method were unacceptably noisy (see the fourth subplots of Fig. 1 and Fig. 2 for typical examples of the RL estimates). For this reason, it was decided to terminate the RL estimations before the steady-state was reached, at the point when the NMSE value was minimized. (Note that such NMSE-optimal RL estimates are impossible to compute in real-life scenarios, as their computation requires knowing the original images.) The SRL algorithm, on the other hand, produced monotonous convergence in NMSE, producing substantially smaller values of NMSE as compared to the case of RL estimation. For the case of high- and low-count scenarios, some typical reconstruction results are depicted in Fig. 1 and Fig. 2, respectively, from which the superiority of the proposed method can be easily appreciated. It is also interesting to note that, as predicted by the theory, the RL algorithm attempts to arrive at a sparse estimate of the signal of interest, which contradicts the piece-wise smooth nature of our test signals. The SRL method, on the other hand, require the sparsity of the representation coefficients, which appears to be a much more realistic assumption for the case at hand.
A quantitative comparison of the two reconstruction algorithms is presented in Fig. 3, which depicts the values of NMSE produced by the RL and SRL methods as a function of the number of iterations. It is important to emphasize that each value of the NMSE in Fig. 3 is a result of averaging the errors obtained in a series of independent trials, where both the true signals and noises were drawn randomly. Observing Fig. 3, one can see that SRL results in considerably lower values of the NMSE as compared to the case of RL. As well, it converges monotonously to a steady-state solution after a relatively small number of iterations. On the other hand, the convergence of RL is not monotone – the fact that necessitates termination of the algorithm after a predefined number of iterations. Needless to say, determining a proper number of iterations is a data-dependent and somewhat esoteric task in practice.
5.3 Two-dimensional Reconstruction
In the second part of the experimental study, the SRL algorithm is examined using imaging data. As reference methods, the standard RL method as well as the RLTV algorithm (as specified by (23)) are employed to assess the performance of SRL in a comparative way. Note that as neither RL nor RLTV could provide acceptable steady-state results, their iterations needed to be terminated at an earlier point. In particular, the iterations were terminated at the point when the NMSE produced by the algorithms reached a minimum value. It is important to reiterate that such “MSE-optimal” termination rule is not realizable in real-life scenarios where the original images are unknown. Hence, the RL and RLTV estimates demonstrated below represent the best possible result which could be achieved using these reconstruction methods.
Implementation of the SRL algorithm relies on the availability of a dictionary of positive-valued atoms . In the 2-D experiments, we explore two possible approaches to the definition of such a dictionary, viz. unsupervised and supervised. In the first case, the dictionary is composed of shift-invariant subsets generated by cubic B-splines , and therefore it does not depend on the nature of imagery data at hand. In the second case, the dictionary is learned from a training set that is supposed to represent the family of objects which the measured image is likely to belong to.
5.3.1 Unsupervised dictionary
For the sake of reproducibility of the results of this paper, we provide some necessary details on the construction of the cubic spline dictionary. To this end, let denote a vector of length , whose elements are all equal to 1. Also, let be the vector of integral values of the cubic B-spline and be the vector of integral values of the cubic B-spline scaled by a dyadic factor of , with . The vectors have the length of points and they can be computed recursively according to 
where denotes the operation of convolution. As prescribed by common practice, the vectors can be subsequently normalized to obey .
The 1-D kernels can be extended into isotropic and separable 2-D (discrete) splines by means of the tensor product, i.e. , . Consequently, the true image is assumed to belong to the space of all integer translations of the discrete kernels . Under periodic boundary conditions, there are a total of representation coefficients for each . Let the array of all these coefficients be denoted by . Then, the operator can be shown to be given by
while its adjoint is given by
with standing for circular convolution111For some examples of related MATLAB codes and their use visit www.ece.uwaterloo.ca/õlegm/research.html.. According to (29), an image produces a total of spline coefficients, and hence the overall complexity of applying and is comparable to that of a stationary wavelet transform.
5.3.2 Supervised dictionary
In the supervised case, the dictionary is designed adaptively based on a set of training examples that represent the family of images, to which is expected to belong. Consequently, the training procedure aims to find a set of positive-valued representation functions, using which one can represent the training images with a minimum possible number of representation coefficients. In this work, the training was performed by means of the the non-negative kSVD (NN-kSVD) algorithm222A software for the NN-kSVD algorithm can be obtained from http://www.cs.technion.ac.il/ẽlad/software/. proposed in . In compliance with the Sparseland model of , the training was applied to overlapping segments of a total of 11 MRI scans of the brain, which did not include the scans used in reconstruction. Each of these segments was represented by a linear combination of a total of 512 atoms, thereby resulting in the over-completeness factor of 2 per segment. For the formal definition of operators and the reader is kindly referreded to .
5.3.3 Reconstruction of MRI scans
The blurring kernel used in our next experiment was defined to be , with , while the regularization parameter was set to be equal to 0.1 in both unsupervised and supervised cases. A typical example of the original image used in this study along with its blurred and noise-contaminated versions are shown in the leftmost, middle, and the rightmost subplots of Fig. 4, respectively. The “measured” MRI scans have been obtained by contaminating the blurred image with Poisson noise giving rise to an SNR of approximately 15 dB.
Some typical reconstruction results produced by the proposed and reference methods are shown in Fig. 5. One can see that the SRL method is capable of better recovering the details of the true image as compared to the alternative approaches. It can also be seen that the SRL estimates possess higher resolution and contrast as compared to the alternative approaches. These conclusions are further supported by the quantitative measures of Table 1, which compares the NMSE and SSIM indices of different estimates. As shown by the table, the supervised SRL method produces the lowest NMSE and the largest SSIM index among all the methods under comparison.
|RL (MSE-optimal)||RLTV (MSE-optimal)||SRL (Splines)||SRL (NN-kSVD)|
6 Discussion and Conclusions
In the current paper, a new regularized version of the RL algorithm has been presented, which is advantageous in a number of ways. First, the proposed method is general in its formulation. The latter allows applying the same reconstruction procedure to a number of different settings, such as image de-noising or image enhancement through deconvolution. Attuned to deal with Poissonian noises, the method can therefore be applied for the reconstruction of imagery data produced by many important image modalities, including optical, microscopic, turbulent, and nuclear imaging, just to name a few. Moreover, whilst many alternative reconstruction methods take advantage of certain simplifying assumptions about the noise nature, the proposed technique is optimized to deal with the realistic noise model at hand.
Another important advantage of the proposed method consists in the generality of the prior model used to describe the nature of true images. Specifically, the images have been assumed to admit a sparse representation in the domain of a properly chosen linear transform. Note that the above a priori model is nowadays considered to be a standard model in numerous applications of the theory of sparse representations. This is what allows the SRL algorithm to recover piecewise smooth images – the task impossible to accomplish by means of the standard RL procedure for its property to favour spatially sparse reconstructions.
Yet another critical advantage of the proposed SRL algorithm is in its algorithmic structure, which allows solving non-smooth optimization problems at the computational cost of a steepest descent procedure. Consequently, the computational load required by the proposed method is relatively small, which allows the method to be applied for the solution of large-scale problems and/or for processing of large data sets.
-  L. A. Shepp and Y. Vardi, “Maximum likelihood reconstruction for emission tomography,” IEEE Trans. Med. Imag., vol. 1, no. 2, pp. 113–122, Oct. 1982.
-  S.-J. Lee, A. Rangarajan, and G. Gindi, “Bayesian image reconstruction in SPECT using higher order mechanical models as priors,” IEEE Trans. Med. Imag., vol. 14, no. 4, pp. 669–680, Dec. 1995.
-  M. Yavuz and J. A. Fessler, “Statistical image reconstruction methods for randomsprecorrected PET scans,” Med. Im. Anal., vol. 2, no. 4, pp. 369–378, 1998.
-  H. H. Bauschke, D. Noll, A. Celler, and J. M. Borwein, “An EM algorithm for dynamic SPECT,” IEEE Trans. Med. Imag., vol. 18, no. 3, pp. 252–261, Mar. 1999.
-  R. A. Boie and I. J. Cox, “An analysis of camera noise,” IEEE Trans. Pattern Anal. Machine Intell., vol. 14, no. 6, pp. 671–674, June 1992.
-  G. E. Healey and R. Kondepudy, “Radiometric CCD camera calibration and noise estimation,” IEEE Trans. Pattern Anal. Machine Intell., vol. 16, no. 3, pp. 267–276, Mar. 1994.
-  T. J. Holmes, “Blind deconvolution of quantum-limited incoherent imagery: Maximum-likelihood approach,” J. Opt. Soc. Am. A, vol. 9, no. 7, pp. 1052–1061, July 1992.
-  H. Bradt, Astronomy methods: A physical approach to astronomical observations. Cambridge University Press, 2004.
-  M. C. Roggeman and B. Welsh, Imaging through turbulence. CRC Press, 1996.
-  F. J. Anscombe, “The transformation of Poisson, binomial, and negative binomial data,” Biometrika, vol. 35, pp. 246–254, 1948.
-  P. Fryzlewicz and G. P. Nason, “A Haar-Fisz algorithm for Poisson intensity estimation,” J. Comput. Graph. Statist., vol. 13, pp. 621–638, 2004.
-  M. Jansen, “Multiscale Poisson data smoothing,” J. Roy. Statist. Soc. B, vol. 68, no. 1, pp. 27–48, 2006.
-  B. Zhang, J. M. Fadili, and J.-L. Starck, “Wavelets, ridgelets, curvelets for Poisson,” IEEE Trans. Image Processing, vol. 17, no. 7, pp. 1093–1108, July 2008.
-  E. D. Kolaczyk, “Nonparametric estimation of intensity maps using Haar wavelets and Poisson noise characteristics,” Astrophys. J., vol. 534, pp. 490–505, 2000.
-  B. Zhang, J. M. Fadili, and J.-L. Starck, “Fast Poison noise removal by biorthogonal Haar domain hypothesis testing,” Statist. Method., no. 5, pp. 387–396, 2008.
-  E. D. Kolaczyk, “Bayesian Multiscale Models for Poisson Process,” J. Amer. Statist. Assoc., vol. 94, no. 447, pp. 920–933, Sept. 1999.
-  R. D. Nowak and R. G. Baraniuk, “Wavelet-domain filtering for photon imaging systems,” IEEE Trans. Image Processing, vol. 8, no. 5, pp. 666–678, May 1999.
-  K. E. Timmermann and R. D. Nowak, “Multiscale modeling and estimation of Poisson processes with application to photo-limited imaging,” IEEE Trans. Inform. Theory, vol. 45, no. 3, pp. 846–862, Apr. 1999.
-  R. D. Nowak and E. D. Kolaczyk, “A statistical multiscale framework for Poisson inverse problems,” IEEE Trans. Inform. Theory, vol. 46, no. 5, pp. 1811–1825, Aug. 2000.
-  S. Sardy, A. Antoniadis, and P. Tseng, “Automatic smoothing with wavelets for a wide class of distributions,” J. Comput. Graph. Statist., vol. 13, no. 2, pp. 399–421, June 2004.
-  A. Antoniadis and T. Sapatinas, “Wavelet shrinkage for natural exponential families with quadratic variance function,” Biometrica, vol. 88, no. 3, pp. 805–820, 2001.
-  W. H. Richardson, “Bayesian-based iterative method of image restoration,” J. Opt. Soc. Am. A, vol. 62, no. 1, pp. 55–59, 1972.
-  L. B. Lucy, “An iterative technique for the rectification of observed distributions,” Astron. J., vol. 79, no. 6, pp. 745–754, 1974.
-  R. Molina, J. Mateos, and J. Abad, “Prior models and the Richardson-Lucy restoration method,” Restoration HST Images Spectra II, vol. 52, no. 118, pp. 118–122, 1994.
-  G. van Kempen anf L.van Vliet, “The influence of the regularization parameter and the first estimate on the performance of Tikhonov regularized non-linear image restoration algorithms.” Journal of Microscopy, vol. 198, no. 1, pp. 63–75, Apr. 2000.
-  N. Dey, L. Blanc-Feraud, C. Zimmer, P. Roux, Z. Kam, J.-C. Olivo-Marin, and J. Zerubia, “Richardson-Lucy algorithm with total variation regularization for 3-D confocal microscope deconvolution,” Miscosc. Res. Techniq., vol. 69, no. 4, pp. 260–266, Apr. 2006.
-  L. Rudin, S. Osher, and E. Fatemi, “Nonlinear total variation based noise removal algorithms,” Physica D: Nonlinear phenomena, vol. 60, no. 1-4, pp. 259–268, 1992.
-  D. L. Donoho, “De-noising by soft-thresholding,” IEEE Trans. Inform. Theory, vol. 41, no. 3, pp. 613–627, May 1995.
-  S. S. Chen, D. L. Donoho, and M. A. Saunders, “Atomic decomposition by Basis Pursuit,” SIAM Review, vol. 43, no. 1, pp. 129–159, Mar. 2001.
-  M. Elad and M. Aharon, “Image denoising via sparse and redundant representations over learned dictionaries,” IEEE Trans. Image Processing, vol. 15, no. 12, pp. 3736–3745, Dec. 2006.
-  M. Elad, B. Matalon, J. Shtok, and M. Zibulevsky, “A wide-angle view at iterated shrinkage algorithms,” in Proceedings of SPIE (Wavelet XII), San-Diego CA, USA, 2007.
-  Z. Wang, A. C. Bovik, H. R. Sheikh, and E. P. Simoncelli, “Image quality assessment: From error visibility to structural similarity,” IEEE Trans. Image Processing, vol. 13, no. 4, pp. 600–612, Apr. 2004.
-  M. Elad, “Why simple shrinkage is still relevant for redundant representations?” IEEE Trans. Inform. Theory, vol. 52, no. 12, pp. 5559–5569, Dec. 2006.
-  D. Donoho, “Compressed sensing,” IEEE Trans. Inform. Theory, vol. 52, no. 4, pp. 1289–1306, Apr. 2006.
-  E. Candes and T. Tao, “Near optimal signal recovery from random projections: Universal encoding strategies,” IEEE Trans. Inform. Theory, vol. 52, no. 12, pp. 5406–5425, Dec. 2006.
-  S. Boyd and L. Vandenberghe, Convex optimization. Cambridge University Press, 2004.
-  P. Bofill and M. Zibulevsky, “Underdetermined blind source separation using sparse representations,” Signal Processing, vol. 81, no. 11, pp. 2353–2362, Nov. 2001.
-  P. Georgiev, F. Theis, and A. Cichocki, “Sparse component analysis and blind source separation of underdetermined mixtures,” IEEE Trans. Neural Networks, vol. 16, no. 4, pp. 992–996, July 2005.
-  I. Daubechies, M. Defrise, and C. DeMol, “An iterative thresholding algorithm for linear inverse problems with a sparsity constraint,” arXiv:math/0307152v2, 2003.
-  M. A. T. Figueiredo and R. D. Nowak, “An EM algorithm for wavelet-based image restoration,” IEEE Trans. Med. Imag., vol. 12, no. 8, pp. 906–916, Aug. 2003.
-  J. L. Starck, M. Elad, and D. Donoho, “Redundant multiscale transforms and their application for morphological component analysis,” Adv. Imag. Electron. Phys., vol. 132, 2004.
-  O. Michailovich and D. Adam, “A high-resolution technique for ultrasound harmonic imaging using sparse representations in Gabor frames,” IEEE Trans. Med. Imag., vol. 21, no. 12, pp. 1490–1503, Dec. 2002.
-  M. Aharon, M. Elad, and A. M. Bruckstein, “K-SVD and its non-negative variant for dictionary design,” in Proceedings of SPIE conference wavelets, vol. 5914, July 2005.
-  M. Unser, “Splines: A perfect fit for signal and image processing,” IEEE Signal Process Mag., vol. 16, no. 6, pp. 22–38, Nov. 1999.