PaperlayerView exact revision on arXiv

Source-linked AI summary

Restoration of Poissonian Images Using Alternating Direction Optimization

Mário A. T. Figueiredo, José M. Bioucas-Dias

arXiv:1001.2244v2math.OCmath.NA

TL;DR

Poissonian image deconvolution is difficult because its likelihood is non-quadratic and non-separable, while regularization is nonsmooth and constrained by non-negativity. The paper applies an ADMM-based method, PIDAL, to TV and frame-based formulations, deriving conditions for solution existence, uniqueness, and convergence. The resulting algorithms are reported to outperform state-of-the-art alternatives in speed and restoration accuracy.

  • Problem

    Poissonian restoration combines a non-quadratic, non-separable likelihood with nonsmooth regularization and a non-negativity constraint.

  • Method

    PIDAL applies an instance of ADMM to Poissonian deconvolution problems with total-variation and frame-based regularization.

  • Results

    The resulting algorithms are reported to outperform state-of-the-art alternatives in speed and restoration accuracy.

  • Takeaways & Limitations

    The paper provides sufficient conditions for solution existence, uniqueness, and convergence across total-variation and frame-based formulations.

Abstract

from arXiv · show

Much research has been devoted to the problem of restoring Poissonian images, namely for medical and astronomical applications. However, the restoration of these images using state-of-the-art regularizers (such as those based on multiscale representations or total variation) is still an active research area, since the associated optimization problems are quite challenging. In this paper, we propose an approach to deconvolving Poissonian images, which is based on an alternating direction optimization method. The standard regularization (or maximum a posteriori) restoration criterion, which combines the Poisson log-likelihood with a (non-smooth) convex regularizer (log-prior), leads to hard optimization problems: the log-likelihood is non-quadratic and non-separable, the regularizer is non-smooth, and there is a non-negativity constraint. Using standard convex analysis tools, we present sufficient conditions for existence and uniqueness of solutions of these optimization problems, for several types of regularizers: total-variation, frame-based analysis, and frame-based synthesis. We attack these problems with an instance of the alternating direction method of multipliers (ADMM), which belongs to the family of augmented Lagrangian algorithms. We study sufficient conditions for convergence and show that these are satisfied, either under total-variation or frame-based (analysis and synthesis) regularization. The resulting algorithms are shown to outperform alternative state-of-the-art methods, both in terms of speed and restoration accuracy.

I. INTRODUCTION

Poissonian image deconvolution combines a difficult convex optimization problem with observation models that differ from Gaussian-noise settings. The paper develops an ADMM-based approach and establishes solution and convergence conditions for several regularizers.

  • I. INTRODUCTION: Poissonian image models are relevant to astronomical, biomedical, and photographic imaging, while algorithms designed for Gaussian noise cannot be directly applied.The observation operator may model convolution, tomography, or denoising through a linear matrix K.
  • I. INTRODUCTION: Poissonian deconvolution uses a convex objective with a smooth but non-quadratic data term, a nonsmooth regularizer, and non-negativity constraints.Its high dimensionality makes direct use of off-the-shelf optimization algorithms difficult.
  • I. INTRODUCTION: Poisson log-likelihood is generally non-separable and lacks a Lipschitz-continuous gradient, limiting guaranteed use of forward-backward methods.When applied despite these difficulties, such methods are slow, especially for severely ill-conditioned observation operators.
  • I. INTRODUCTION: The paper proposes PIDAL, an image-deconvolution algorithm based on an instance of the alternating direction method of multipliers.The method belongs to the augmented Lagrangian family and uses a splitting strategy distinct from split-Bregman.
  • I. INTRODUCTION: The study covers total variation, frame-based analysis, and frame-based synthesis regularization, with sufficient conditions for existence and uniqueness of solutions.It also develops convergence conditions and evaluates the resulting algorithms against state-of-the-art alternatives.
  • I. INTRODUCTION: The negative log-likelihood is constructed from component functions ξ(z, y) that are convex and coercive, and strictly convex when y > 0.For L and L ◦ K, strict convexity and coercivity additionally depend on positive observations and injectivity of K.

B. Regularization Criteria: Analysis and Synthesis Formulations

The paper formulates Poissonian restoration with MAP objectives combining likelihood, regularization, and non-negativity constraints. It analyzes total variation and frame-based priors, deriving conditions under which minimizers exist or are unique.

  • B. Regularization Criteria: Analysis and Synthesis Formulations: MAP restoration minimizes a Poisson log-likelihood plus a regularizer under a positivity constraint.The paper introduces three corresponding objective constructions.
  • 1) Total Variation Regularization:: Isotropic discrete total variation regularization penalizes first-order horizontal and vertical image differences.Its optimization problem includes the indicator of the non-negative orthant.
  • 1) Total Variation Regularization:: For TV regularization, the objective is proper, lower semicontinuous, and convex.These properties follow from the corresponding properties of the likelihood, TV term, and non-negativity indicator.
  • 1) Total Variation Regularization:: TV minimizers exist when the constant-image null space intersects the null space of K only at zero.Existence also holds when K is injective, or when K is nonnegative with at least one strictly positive element.
  • 1) Total Variation Regularization:: TV has a unique solution when K is injective and every observed count y_i is nonzero.Under these conditions, the objective is coercive and strictly convex.
  • 2) Frame Analysis Regularization:: Frame-based analysis regularization penalizes representation coefficients Px, typically with an ℓ1 norm.Here P is the analysis operator associated with a wavelet basis or tight frame.

2) Frame Analysis Regularization:

Frame analysis regularization applies an analysis operator before penalizing coefficients, while enforcing non-negativity on the image estimate. The resulting criterion always has a minimizer and is unique under injectivity conditions.

  • Frame analysis represents regularization through an analysis operator, with non-negativity imposed on the estimated image.
  • LFA is proper, lower semi-continuous, convex, and coercive, so it has a minimizer.
  • If K is injective and every observation count yi is nonzero, LFA is strictly convex with a unique minimizer.
  • Frame synthesis instead estimates frame coefficients s in x = W s, while applying non-negativity to W s rather than to the coefficients.
  • For frame synthesis, LFS has a minimizer, and it is unique when K W is injective and every yi is nonzero.

III. THE ALTERNATING DIRECTION METHOD OF MULTIPLIERS (ADMM)

ADMM solves a constrained splitting of convex objectives by alternating minimizations and dual updates. A multi-function mapping preserves separability and convergence under full-column-rank and summable-error conditions.

  • The standard ADMM: ADMM is introduced for minimizing f1(z) + f2(Gz), with f1 and f2 closed, proper, convex functions.
  • The standard ADMM: Under full column rank of G and existence of a solution, the ADMM primal sequence converges to a solution.
  • The standard ADMM: Exact minimization is unnecessary when the corresponding error sequences are absolutely summable, because convergence is preserved.
  • A variant of ADMM: The proposed multi-function mapping represents a sum of J convex terms within the two-function ADMM form.
  • A variant of ADMM: With f1 = 0, the z-update becomes a unique quadratic minimization when G has full column rank.
  • A variant of ADMM: The chosen mapping decouples the u-updates into J independent minimizations, including Moreau proximity operators for the individual functions.
  • A variant of ADMM: For the resulting ADMM, convergence is guaranteed when G has full column rank and all constituent functions are closed, proper, and convex.

A. Applying ADMM

Applying the multi-function ADMM mapping to total variation yields PIDAL-TV, whose updates combine quadratic, separable, denoising, and non-negativity steps. Under periodic convolution, its computational cost scales as O(n log n).

  • Applying ADMM: The framework is applied to the total-variation-based Poisson deconvolution criterion.
  • Applying ADMM: The resulting algorithm, PIDAL-TV, is an ADMM method for Poisson image deconvolution with total variation regularization.
  • Implementation aspects: PIDAL-TV uses G = [K^T I I]^T, which has full column rank; periodic convolution enables its inversion through FFT operations.
  • Implementation aspects: The TV update is a Moreau proximity operator corresponding to TV denoising, while the non-negativity update projects onto the first orthant.
  • Computational cost: Under periodic convolution, FFT-based products cost O(n log n), while the remaining listed updates have O(n) cost.

C. Convergence of PIDAL-TV: Exact TV

PIDAL-TV converges to a minimizer when its TV subproblem is solved exactly, and approximate TV updates can retain convergence when their errors are summable. Warm-starting Chambolle’s algorithm supports this condition with few inner iterations.

  • Exact TV: If the minimizations in lines 9, 11, and 13 are solved exactly, PIDAL-TV converges to a minimizer of (10), provided one exists.
  • Approximate TV: The TV subproblem has no closed-form solution, so PIDAL-TV uses Chambolle’s iterative denoising algorithm.
  • Approximate TV: Approximate TV minimization does not compromise ADMM convergence when its error sequence is summable.
  • Approximate TV: Warm-starting Chambolle’s internal variables from the previous PIDAL-TV iteration is used to obtain summable errors.
  • Approximate TV: The fitted error tails support summability, and experiments therefore use t = 5 Chambolle iterations per PIDAL-TV call.
  • Approximate TV: For t = 5 and t = 20 Chambolle iterations, the error sequences ρ(a)_k approach zero, whereas ρ(b)_k sequences are not decreasing.

V. POISSONIAN IMAGE RECONSTRUCTION WITH FRAME-BASED REGULARIZATION

The paper develops PIDAL algorithms for frame-based analysis and synthesis regularization, with convergence guarantees and O(n log n) computational scaling for the analysis variant.

  • Frame-based analysis: The analysis formulation uses H^(1) = K, H^(2) = P, and H^(3) = I within the ADMM construction.
  • Frame-based analysis: PIDAL-FA applies ADMM to the frame-based analysis criterion and uses a Parseval-frame assumption to simplify its matrix inversion.For convolution operators, the inverse can be computed with FFTs at O(n log n) cost.
  • Frame-based analysis: The ℓ1 minimization step in PIDAL-FA is implemented by the Moreau proximity operator, corresponding to soft-thresholding.
  • Frame-based analysis: O(n log n) is the overall computational cost of PIDAL-FA.Its O(n log n) steps include the main matrix operations, while several other lines have O(n) cost.
  • Convergence: PIDAL-FA converges to a minimizer of the frame-based analysis objective, provided one exists.

B. Synthesis Criterion

PIDAL-FS is the ADMM algorithm for frame-based synthesis regularization. Its matrix inversion is accelerated using Parseval-frame structure and the Sherman–Morrison–Woodbury identity, while convergence is guaranteed when a minimizer exists.

  • Synthesis algorithm: The synthesis construction uses H^(1) = KW, H^(2) = I, and H^(3) = W.
  • Synthesis algorithm: PIDAL-FS applies ADMM to the frame-based synthesis objective and is presented as the resulting synthesis algorithm.
  • Efficient implementation: For a 1-tight frame, PIDAL-FS uses the Sherman–Morrison–Woodbury identity to avoid directly computing the difficult matrix inverse.The difficulty arises because KW is not block circulant even when K represents a periodic convolution.
  • Efficient implementation: The relevant inversions have O(n) cost, while FFT-based products by the frame-related operators contribute O(n log n) cost.
  • Computational cost: O(n log n) is the overall computational cost of PIDAL-FS.The remaining listed operations have O(n) cost.
  • Convergence: PIDAL-FS converges to a minimizer of the synthesis objective, provided one exists.The convergence argument uses exact closed-form minimizations and full column rank of the associated matrix.

VI. EXPERIMENTS

The experiments compare PIDAL methods with other state-of-the-art approaches under MATLAB implementations and tune regularization for maximum ISNR improvement. The ADMM parameter is chosen heuristically for satisfactory speed.

  • Experimental setup: PIDAL is compared with state-of-the-art methods proposed in references,, and.All algorithms are implemented in MATLAB on a 3.0GHz Intel Core2Extreme PC with 4GB RAM under Windows Vista.
  • Experimental setup: The regularization parameter τ is adjusted to achieve the highest improvement in signal-to-noise ratio.The experiments define ISNR as 10 log10 of the relevant signal-to-noise ratio improvement.
  • Parameter choice: ADMM converges for any µ > 0, but µ affects algorithm speed.The experiments use the rule of thumb µ = 60τ/M, where M is the maximum intensity of the original image.
  • Parameter choice: Changing µ by one order of magnitude below or above the selected rule does not significantly change the reported results.

A. Comparison with [37]

The experiments compare PIDAL-TV, PIDAL-FA, and PIDAL-FS with competing Poisson-image restoration methods across several image settings. PIDAL-TV and PIDAL-FA generally provide strong restoration accuracy and computational efficiency, while PIDAL-FS is consistently weaker in the reported comparisons.

  • Comparison with [37]: PIDAL-TV is compared with PIDSplit and PIDSplit+ on a blurred, noisy Cameraman image using matched iteration budgets.PIDAL-TV uses five Chambolle iterations per outer iteration and is run for 430 iterations to approximate the competing methods’ 2150 iterations.
  • Frame-based comparisons: PIDAL-FA slightly outperforms PIDAL-TV in ISNR, has similar MAE, and reaches its best estimate faster in the reported example.PIDAL-FS is slightly worse in both ISNR and MAE and takes longer than PIDAL-FA; this pattern recurs in the later experiments.
  • Convergence: PIDAL-TV reaches convergence after about 160 iterations, corresponding to 25 seconds, in the high-SNR Cameraman experiment.The experiment tracks MAE and ISNR over iterations and elapsed time.
  • Experimental design: The experiments evaluate PIDAL-TV and PIDAL-FA against state-of-the-art methods using ISNR, MAE, and running time.The comparisons include DFS and RL-MRS, with image intensities scaled to M ∈ {5, 30, 100, 255}.
  • Overall findings: The paper’s overall experiments report state-of-the-art performance in both speed and restoration accuracy.The authors frame the contribution as algorithms for solving variational Poisson-restoration optimization problems rather than as a new restoration criterion.

APPENDIX A: CONVEX ANALYSIS

The appendix reviews convex-analysis concepts used to establish optimization results. It defines convexity and related function properties, then states existence and uniqueness principles for minimizers.

  • Convexity: A function is convex when its value at a convex combination does not exceed the corresponding combination of function values.Strict convexity replaces the inequality with a strict inequality for distinct points in the domain and interior combination weights.
  • Function properties: A function is proper when it is not everywhere equal to +∞ and never equals −∞.This property excludes degenerate extended-real-valued functions.
  • Function properties: A function is coercive when its value tends to +∞ as the norm of its argument tends to infinity.Lower semicontinuity is defined through local behavior on norm balls.
  • Minimization theorems: A proper, lower-semicontinuous, coercive, convex function has a nonempty set of minimizers.This theorem supplies an existence condition for optimization solutions.
  • Minimization theorems: A strictly convex function has at most one minimizer.Strict convexity therefore provides a uniqueness condition when a minimizer exists.
Loading 1001.2244v2…