PaperlayerView exact revision on arXiv

Source-linked AI summary

A proximal iteration for deconvolving Poisson noisy images using sparse representations

François-Xavier Dupé, Jalal Fadili, Jean Luc Starck

arXiv:0803.2623v2stat.APmath.OCmath.ST

TL;DR

The paper tackles image deconvolution when observations are blurred and contaminated by Poisson noise, especially when the target image is sparsely represented. It uses a variance-stabilizing transform and sparse convex optimization solved by a fast proximal iteration. Experiments report benefits from modeling Poisson statistics and image morphology, including encouraging results on real confocal microscopy images.

  • Problem

    The problem is to restore an image x from an observed count image y when deconvolution data contain Poisson noise and sparsity-based methods have mainly addressed Gaussian noise.

  • Method

    The method applies the Anscombe variance stabilizing transform and minimizes a convex noise-aware functional with sparse representation penalties and positivity constraints using a fast proximal iterative algorithm.

  • Results

    The experiments report striking benefits from accounting for Poisson statistics and image morphology, with encouraging preliminary results on real confocal microscopy images.

  • Takeaways & Limitations

    Sparse-domain regularization may be tractable for Poisson-noise deconvolution applications such as astronomy and microscopy.

Abstract

from arXiv · show

We propose an image deconvolution algorithm when the data is contaminated by Poisson noise. The image to restore is assumed to be sparsely represented in a dictionary of waveforms such as the wavelet or curvelet transforms. Our key contributions are: First, we handle the Poisson noise properly by using the Anscombe variance stabilizing transform leading to a {\it non-linear} degradation equation with additive Gaussian noise. Second, the deconvolution problem is formulated as the minimization of a convex functional with a data-fidelity term reflecting the noise properties, and a non-smooth sparsity-promoting penalties over the image representation coefficients (e.g. $\ell_1$-norm). Third, a fast iterative backward-forward splitting algorithm is proposed to solve the minimization problem. We derive existence and uniqueness conditions of the solution, and establish convergence of the iterative algorithm. Finally, a GCV-based model selection procedure is proposed to objectively select the regularization parameter. Experimental results are carried out to show the striking benefits gained from taking into account the Poisson statistics of the noise. These results also suggest that using sparse-domain regularization may be tractable in many deconvolution applications with Poisson noise such as astronomy and microscopy.

I. INTRODUCTION

The paper addresses image deconvolution under Poisson noise using sparse representations and a noise-aware optimization framework. It combines the Anscombe transform, convex sparsity-regularized minimization, and a fast proximal algorithm, with experiments showing benefits from modeling Poisson statistics and image morphology.

  • Sparsity-based deconvolution methods had mainly focused on Gaussian noise, leaving Poisson-noise deconvolution comparatively under-addressed.
  • The paper proposes deconvolution for blurred images contaminated by Poisson noise, targeting applications including biomedical imaging and astronomy.
  • The Anscombe variance stabilizing transform converts the Poisson observation model into a nonlinear degradation equation with additive Gaussian noise.
  • The restoration minimizes a convex functional combining noise-reflecting nonlinear data fidelity, sparsity-promoting coefficient penalties, and positivity constraints.The representation may use orthogonal wavelets, translation-invariant wavelets, curvelets, or combined wavelet–curvelet systems.
  • A fast proximal iterative algorithm solves the minimization problem while supporting overcomplete sparse multiscale representations of image morphology.
  • Experiments on simulated and real images report striking benefits from accounting for Poisson noise and morphological structures through overcomplete sparse multiscale transforms.

A. Relation to prior work

The paper addresses Poisson deconvolution by combining variance stabilization with an optimization formulation that preserves the noise statistics and supports proximal splitting.

  • Motivation: Gaussian-noise deconvolution methods are inadequate because Poisson noise becomes approximately Gaussian only at large intensities and has intensity-dependent variance.The paper therefore treats the Poisson statistics explicitly rather than relying on a conventional Gaussian model.
  • Optimization formulation: A Poisson likelihood-based fidelity term reflects the conditional distribution of observed counts given the original image.This corresponds to an anti-log-likelihood score derived from the acquisition model.
  • Algorithmic motivation: Direct MAP optimization is incompatible with forward-backward splitting because the resulting data-fidelity gradient is not Lipschitz-continuous.The paper motivates the transformed formulation partly because Lipschitz continuity is needed for uniquely defined iterates and convergence.
  • Anscombe modeling: The Anscombe variance-stabilizing transform produces additive unit-variance Gaussian noise but makes the degradation equation nonlinear.The nonlinearity arises from the transform's square-root relationship.
  • Contribution: The proposed work formulates and solves the nonlinear Poisson problem with an optimization and fixed-point algorithm, reporting benefits from modeling Poisson noise.The paper places this contribution within a broader treatment of sparse representations and monotone-operator splitting.

III. SPARSE IMAGE REPRESENTATION

The image is represented sparsely using an orthogonal basis or tight-frame dictionary, with transform choice guided by the morphology of image content and restoration needs.

  • Sparse representation: Images are modeled as superpositions of elementary atoms and sought in representations containing only a few significant coefficients.Wavelet, curvelet, contourlet, wave-atom, and local-DCT systems are cited as effective sparse representations for different image structures.
  • Dictionary construction: The dictionary is an n × L matrix of unit-norm waveforms, with L ≥ n and L > n indicating a redundant or overcomplete representation.Its analysis transform is T = ΦT.
  • Dictionary construction: The paper restricts the dictionary to an orthobasis or tight frame, for which reconstruction has simple algebraic forms.For a tight frame, ΦΦT = cI and pseudo-inverse synthesis reduces to c^-1Φ.
  • Dictionary design: Morphological diversity combines transforms so each sparsely represents the image component it serves while inefficiently representing other content.This motivates combining complementary systems such as wavelets for isotropic singularities and curvelets for contours.
  • Optimization formulation: The data-fidelity component is derived from the transformed observation model and is convex with a Lipschitz-continuous gradient.Its strict convexity requires an orthobasis with ker(H) = ∅.
  • Optimization formulation: The optimization penalizes representation coefficients to promote sparsity, while an indicator term imposes positivity on the restored image.The positivity constraint reflects the fact that fitted Poisson intensities are positive.

B. Characterization of the solution

The paper characterizes the convex minimization problem, decomposes its nonsmooth terms through proximity operators, and solves it with convergent proximal iterations under stated conditions.

  • Existence: The objective has at least one solution because it is coercive and convex.The paper denotes the solution set by M and establishes nonemptiness.
  • Uniqueness: The solution is unique if Φ is an orthobasis and ker(H) = ∅, or if the penalty ψ is strictly convex.Under the orthobasis condition, the data term is strictly convex; strict convexity of ψ instead makes the regularization term strictly convex.
  • Proximity operators: The nonsmooth term combines positivity and sparsity, so its proximity operator generally lacks an explicit form when Φ ≠ I.The paper replaces direct evaluation with a splitting sequence that separately activates the positivity projection and sparsity proximity operator.
  • Proximity operators: Douglas-Rachford splitting handles the two nonsmooth components, whose proximity operators have closed-form expressions under the stated construction.For an orthobasis, the positivity projection in coefficient space is obtained by transforming to the image domain, projecting, and transforming back.
  • Penalty choice: For the ℓ1 penalty, the coefficient proximity operator is component-wise soft-thresholding with threshold λ.This gives a concrete sparse-regularization instance of the general penalty framework.
  • Convergence: Under bounded positive step sizes, relaxation parameters, and controlled errors, the forward-backward iterates converge to a solution.The error sequences support robustness to numerical errors when computing the gradient and proximity operator.

D. Choice of µ

The section addresses selecting the relaxation parameter and regularization parameter, emphasizing convergence conditions and data-driven GCV selection. The proposed degrees-of-freedom approximation makes GCV readily available despite analytical and computational challenges.

  • Choice of µ: The relaxation parameter µ strongly affects convergence speed, and an Armijo–Goldstein strategy adaptively updates µt.The update uses the largest admissible step satisfying the prescribed condition; τ = 1/2 is typical.
  • Convergence: The iterates αt converge to a minimum of J under the stated conditions on the initial point and relaxation sequence.The theorem provides a sufficient convergence condition based on the Lipschitz constant, which may be pessimistic in some applications.
  • Choice of µ: For tight frames, computing proxf2 is simpler because the Douglas–Rachford sub-iteration is unnecessary, though other transforms incur additional cost.The simplified computation comes at the price of extra checks or computational effort for general transforms.
  • GCV model selection: GCV selects λ by balancing sparsity regularization and deconvolution data fidelity, avoiding cumbersome brute-force testing and visual assessment.The criterion is adaptive and data-driven; AIC and BIC are also identified as possible alternatives.
  • GCV model selection: The effective degrees of freedom is difficult to derive because the observation model and the iterative solution are nonlinear; bootstrap and perturbation methods are computationally prohibitive.The paper therefore conjectures that the support cardinality of α⋆ provides a simple estimator of degrees of freedom, making GCV readily available.
  • GCV model selection: The GCV minimum occurs close to the MAE and MSE minima in experiments, although its selected λ is slightly higher than the corresponding optimum.This approximation performed reasonably well across the reported experiments.

F. Computational complexity and implementation details

The implementation relies on fast implicit transform and convolution operators, with computational cost determined by the dictionary, image size, and iteration counts. Orthobases reduce transform computations, while one algorithmic variant avoids choosing the relaxation parameter.

  • Operators and costs: The main computation applies Φ, ΦT, H, and HT as fast implicit operators rather than constructing them explicitly.Applying H or HT uses two FFTs, costing 2n log n operations for n pixels.
  • Operators and costs: Transform complexity depends on the dictionary: orthogonal wavelets cost O(n), whereas TI-DWT and curvelets cost O(n log n).The convolution operators have their own FFT-based cost, and total complexity also depends on L, NFB, and NDR.
  • Implementation choices: With an orthobasis, Φ is bijective, allowing the optimization to be expressed in image samples and reducing computations in iterations (13) and (14).The orthobasis case requires fewer multiplications by Φ and ΦT.
  • Implementation choices: GCV first selects λ with the forward–backward algorithm, after which iteration (14) is used to avoid user selection of the relaxation parameter µ.The second algorithm is chosen after regularization selection specifically to exempt the user from choosing µ.

A. Simulated data

Simulated experiments across neuron, cell, Cameraman, sky, and mixed-structure images evaluate visual quality, intensity-dependent error, dictionary choice, and execution time. The proposed method generally performs well, especially at low and medium intensities, while competing methods have image- or regime-specific strengths.

  • Neuron and low-intensity results: At maximum intensity 30, NaiveGauss performs poorly because its Gaussian noise model under-regularizes the estimate and leaves signal-dependent Poisson noise uncontrolled.RL-TV deconvolves effectively but produces background artifacts and staircase-like structures, whereas the proposed method restores the neuron spine with a visually pleasant result.
  • Quantitative comparisons: The proposed algorithm performs very well across intensity regimes, especially medium to low, while NaiveGauss is among the worst at low intensity.RL-MRS is effective at low and medium intensities and outperforms the proposed method on the Cell image; RL-TV becomes competitive at high intensity.
  • Cell image: For the simulated cell, the proposed method produces a sharper result than FTITPR, while RL-TV restores small details but introduces dominating staircase artifacts.The visual comparison agrees with the MAE results, with RL-MRS showing the best behavior on this image.
  • Cameraman image: For Cameraman at maximum intensity 30, the proposed method gives the most visually pleasing balance between regularization and deconvolution, despite persistent artifacts.RL-MRS has more artifacts and loses photometry, FTITPR is oversmooth with missing details, and RL-TV and NaiveGauss produce many artifacts.
  • Simulated sky: For the simulated sky at maximum intensity 5000, RL-MRS performs best, while the proposed method preserves faint objects but has a less clean background.At this high intensity, NaiveGauss is comparable to the proposed method on galaxies, while FTITPR loses faint objects.
  • Dictionary influence: Dictionary choice depends on image morphology: TI-DWT favors isotropic or localized structures, curvelets favor lines, and combining both can improve mixed content.For the LinesGaussians image, the combined dictionary visually takes the best of both transforms; the best dictionary varies across the Neuron phantom, Cell, and LinesGaussians images.

B. Real data

The method is applied to a real confocal microscopy image of neurons using a wavelet transform and a modeled optical point-spread function. The restored image cleans the background and recovers neuronal structures, but some information is lost.

  • Real microscopy experiment: On the real 512 × 512 confocal neuron image, the method uses a Gaussian approximation of the microscope’s optical PSF and a wavelet transform.The images are displayed in log-scale for visual rendering.
  • Real microscopy experiment: The restoration cleans the background, makes some structures reappear, restores spines, and reconstructs part of the dendritic tree.These effects are reported from the restored image relative to the observed noisy image.
  • Limitations: Some information can be lost, including tiny holes, and the result may improve with a more accurate PSF model.This identifies PSF accuracy as a practical boundary of the reported real-data result.

C. Reproducible research

The authors provide a freely downloadable toolbox containing Matlab functions, scripts, and datasets for reproducible Poisson-noise image deconvolution experiments.

  • The toolbox is freely available for download from the first author’s webpage.
  • It includes Matlab functions, scripts, and datasets for image deconvolution under Poisson noise.
  • The toolbox requires at least WaveLab 8.02.
  • Its scripts reproduce most of the figures included in the paper.

VI. CONCLUSION

The paper develops a sparsity-based iterative deconvolution method for Poisson-noisy images, with principled optimization, parameter selection, and encouraging experimental performance. The authors identify extensions and unresolved analytical issues for future work.

  • The proposed algorithm combines sparse-domain regularization with fast iterative thresholding for image deconvolution under Poisson noise.The approach is designed for data affected by Poisson noise and uses a fast proximal iterative solver.
  • The Poisson noise is handled using a variance-stabilizing approach, while the optimization problem and iterative algorithm receive a theoretical characterization.
  • A GCV-based procedure is used to select the regularization parameter.
  • The method compares favorably with some state-of-the-art algorithms and produces encouraging preliminary results on real confocal microscopy images.
  • The approach generalizes to other nonlinearities when the corresponding data-fidelity term is convex with a Lipschitz-continuous gradient.
  • The degrees-of-freedom expression lacks a rigorous proof, and extensions to 3D confocal volumes and multivalued images remain under investigation.

APPENDIX

The appendix establishes convexity, solution existence and uniqueness conditions, and convergence properties for the optimization and splitting procedures used by the deconvolution method.

  • The data-fidelity component is convex, with strict convexity under an orthobasis and a trivial null space of H.
  • The gradient computation is straightforward, and its Lipschitz continuity follows from boundedness assumptions on the relevant operators and variables.
  • The objective has a solution because it is coercive; uniqueness follows when the dictionary is an orthobasis with ker(H)=∅ or when ψ is strictly convex.
  • The proximal formulation identifies proxλΨ(α) as the unique minimizer of a quadratic perturbation of λΨ.
  • In the constrained case, prox_f2(α) equals proxλΨ(α) when α belongs to C′.
  • The general nonsmooth case uses Douglas–Rachford splitting, and the iterates converge to a unique point under the stated conditions.
  • The forward-backward convergence result follows by combining the cited theorem with the preceding lemmas and proposition.
Loading 0803.2623v2…