Reprint of Inexact Bregman iteration for deconvolution of superimposed extended and point sources,☆☆

https://doi.org/10.1016/j.cnsns.2014.10.020Get rights and content

Highlights

  • The deconvolution of high contrast images consisting of very bright stars and smooth structures around the stars is studied.

  • To restore the region around the stars, the object can be regarded as superposition of point source and extended source.

  • When the position of the point source is known, we introduce a regularization term only for the diffuse component.

  • We give conditions for the solution of the related variational problem for Poisson data with Tikhonov-like regularization.

  • In presence of an overestimation of the regularization parameter, we solve by the inexact Bregman iteration with SGP method.

Abstract

In this paper we consider the deconvolution of high contrast images consisting of very bright stars (point component) and smooth structures underlying the stars (diffuse component). A typical case is a weak diffuse jet line emission superimposed to a strong stellar continuum. In order to reconstruct the diffuse component, the original object can be regarded as the sum of these two components. When the position of the point sources is known, a regularization term can be introduced for the second component. An approximation of the original object can be obtained by solving a reduced variational problem whose unknowns are the intensities of the stars and the diffuse component. We analyze this problem when the detected image is corrupted by Poisson noise and Tikhonov-like regularization is used, giving conditions for the existence and the uniqueness of the solution. Furthermore, since only an overestimation of the regularization parameter is available, we propose to solve the variational problem by inexact Bregman iteration combined with a Scaled Gradient Projection method (SGP). Numerical simulations show that the images obtained with this approach enable us to reconstruct the original intensity distribution around the point source with satisfactory accuracy.

Introduction

Deblurring images corrupted by Poisson noise has recently been subject of an increasingly amount of works, due to the interest in biological and astronomical imaging (see, for example [1], [2], [3], [4], [5], [6], [7] and references therein).

In the framework of high resolution imaging from ground-based telescopes, the deconvolution and reconstruction process can be particularly challenging for high contrast images, i.e. when one is interested in detecting a weak signal close to a very bright star. This is the case, for instance, of the study of circumstellar structures around Young Stellar Objects (YSO). Previous studies on YSO simulations [8] shows that the reconstruction of the diffuse component by standard deconvolution methods is unsatisfactory in the region close to the emitting star. To overcome this difficulty, the basic idea is to consider the original object x as a superposition of two components, a point source xP and an extended source xEx=xP+xEand to introduce different regularizations for the two components. A short review of the methods based on this approach for Gaussian data (least-squares formulation) is provided by [9]. In the framework of Poisson noise, the same idea was proposed by Lucy in [10] and, recently in [11], [12]. In these last papers, the authors assume that the point source is a single star with known position and they propose to introduce a regularization term only for the extended source. Then, in the Bayesian framework, a maximum-likelihood solution can be obtained by minimizing with respect to xP and xE the following constrained nonlinear programming problemminxP,xE0f0(g;H(xP+xE)+b)+βf1(xE)where gRm is the detected image, HRm×n is the imaging matrix, modeling the distortion due to the acquisition system, and b is a given nonnegative constant background term. The objective function is the combination of two terms: the first one measures the discrepancy from the data g and, for Poisson noise, is expressed by the generalized Kullback–Leibler (KL) divergence (or Csiszár I-divergence) given byf0(g;Hx+b)=i=1mgiloggi(Hx+b)i+(Hx+b)i-giwith the agreement 0log0=0. The second function f1 is a regularization term, that enables to incorporate a priori information about the expected diffuse component of the solution into the reconstruction process. f1 is a convex nonnegative function, weighted by the positive regularization parameter β, which has to be estimated. In [11], a Tikhonov-like regularization is used for the diffuse component and a scaled gradient method is proposed to obtain an approximate solution of (1). This approach depends on the choice of the regularization parameter β. In the case of regularization with Poisson data, this selection is basically an open problem. Recently, few discrepancy principles were proposed and discussed (see for example [7], [2], [13]); further, in [14], [15], a constrained model for seminorm regularization is considered when some knowledge of level noise γ is available; this model allows to compute an estimate of β and a restored image by the minimization of a seminorm, under the constraint that the data fidelity is not greater than γ. But also for this formulation an accurate estimate of γ must be known.

A different approach is based on Bregman regularization. The Bregman approach allows to use an overestimated value of β and, in the framework of inverse problems, this procedure has the typical semi-convergence property of the iterative methods, described, for instance, in [16]. The procedure is proposed by Osher et al. in [17] to recover images corrupted by Gaussian noise and, later, by Brune et al. in [18] for dealing with Poisson data. Each iteration of this scheme requires to exactly compute the minimizer of a problem similar to (1), where f1 is replaced by its Bregman distance at the current iterate. When the data fidelity function is the KL divergence, a closed formula for the minimizer may be unavailable and it could be extremely expensive to obtain exact solutions of these subproblems. On the other hand, in recent years the variational model (1) has been deeply investigated in order to design efficient iterative algorithms specifically tailored for different noise statistics and different regularization terms. Then in [19] an inexact version of the iterative procedure is proposed, where the inexactness in the inner subproblem solution of the Bregman iteration is controlled by a criterion that preserves its convergence and its features in image restoration problems.

In this work, we propose to apply the inexact Bregman iteration for the deconvolution of high contrast images with Poisson noise, where the regularization term depends only on the diffuse component and an overestimation of β is used; for the effectiveness of the procedure, the crucial issue becomes to devise an efficient algorithm for approximately solving the inner subproblems. For Tikhonov regularization, the iterates of the inexact algorithm can be efficiently obtained by the so-called Scaled Gradient Projection (SGP) method [20], but the approach can be extended to other regularization terms as, for example, the edge-preserving regularization [21]. For seminorm regularization, suitable inner solvers can be devised so that the convergence of the inexact Bregman approach is assured (for example [22], [23]).

In order to show the effectiveness of the inexact Bregman iteration, we report the analysis of the numerical reconstructions obtained on a simulation of observations of a YSO jet [11]. We consider different star-jet intensity contrast cases by varying the magnitude of the star at each simulation. The obtained results appear very promising for the application to real cases.

The paper is organized as follows. In Section 2, we introduce the features of the problem (1), (2)) and we derive a reduced formulation. We give the condition for the existence and the uniqueness of the solution for Tikhonov-like regularization term and we discuss the implementation of SGP method for its numerical solution. In Section 3, we recall the exact and inexact Bregman method ant its convergence properties and we discuss the application of the scheme for the problem (1), (2). Finally in Section 4 we describe the results obtained by the inexact Bregman scheme on test-problems generated as described in [11] and kindly provided by the authors of this last paper.

Section snippets

Statement of the problem

A common problem in astrophysical imaging is to reconstruct two maps of a given region: one consisting of stars (point sources) and the other representing the smooth structures underlying the stars (extended sources). In this case, the object x to be reconstructed can be considered as the sum of two nonnegative components, corresponding to the point sources (xP) and to the extended sources (xE). Following several recent papers [24], [9], [25], the two components can be reconstructed using

Exact and inexact Bregman procedure for deblurring of Poisson data

In this section we resume the features of the exact and inexact Bregman procedure in the image restoration framework.

First of all, we recall that the Bregman distance of a proper convex function F:RnR between x and y is defined as followsDpF(x,y)=F(x)-F(y)-p,x-ywhere p is a subgradient of F at y and ·,· denotes the canonical inner product of two vectors of Rn.

Numerical experiments

This section is devoted to numerically evaluate the effectiveness of the inexact procedure based on Bregman iteration for recovering high contrast images corrupted by Poisson noise. The numerical experiments described in this section have been performed in MATLAB environment (R2012b), on a PC equipped with an Intel Core i7-3517U processor 1.9 GHz, 8 GB RAM.

We consider four test-problems, provided by the authors of [12], that are simulated Large Binocular Telescope (LBT) infrared narrow-band

Conclusions

In this paper we consider the deconvolution problem of high contrast images corrupted by Poisson noise. These images consist of very bright stars and other smooth structures underlying the stars. In order to investigate the regions around the stars, the original object is written as sum of two components, corresponding to the point sources xP and to the extended sources xE. When the point sources have known positions, the unknowns can be reduced to the intensities of the stars and the pixels of

Acknowledgments

The authors thank Prof. M. Bertero for helpful suggestions and interesting discussions about the variational model and the numerical results; the authors are grateful to Dr. A. La Camera and Prof. P. Boccacci for providing the test-problems used in the simulations of the last section.

References (35)

  • S. Setzer et al.

    Deblurring Poissonian images by split Bregman techniques

    J Vis Commun Image R

    (2010)
  • L.M. Bregman

    The relaxation method of finding the common points of convex sets and its applications to the solution of problems in convex optimization

    USSR Comput Math Math Phys

    (1967)
  • M. Bertero et al.

    Image deblurring with Poisson data: from cells to galaxies

    Inverse Prob

    (2009)
  • J.M. Bardsley et al.

    Regularization parameter selection methods for ill-posed Poisson maximum likelihood estimation

    Inverse Prob

    (2009)
  • R. Zanella et al.

    Efficient gradient projection methods for edge-preserving removal of Poisson noise

    Inverse Prob

    (2009)
  • C. Brune et al.

    Scale space and variationale methods in computer vision

    (2009)
  • M.A.T. Figueiredo et al.

    Restoration of Poissonian images using alternating direction optimization

    IEEE Trans Image Process

    (2010)
  • M. Bertero et al.

    A discrepancy principle for Poisson data

    Inverse Prob

    (2010)
  • P. Ciliegi et al.

    Analysis of LBT LINC-NIRVANA simulated images of galaxies and young stellar objects

  • J.F. Giovanelli et al.

    Positive deconvolution for superimposed extended source and point sources

    Astron Astrophys

    (2005)
  • L.B. Lucy

    An iterative technique for the rectification of observed distributions

    Astron J

    (1974)
  • A.L. Camera et al.

    Image reconstruction for observations with an high dynamic range: LINC-NIRVANA simulations of a stellar jet

  • A.L. Camera et al.

    Reconstruction of high dynamic range images: Simulations of LBT observations of a stellar jet, a pathfinder study for future AO-assisted giant telescopes

    Publ Astron Soc Pac

    (2014)
  • A. Staglianò et al.

    Analysis of an approximate model for Poisson data reconstruction and a related discrepancy principle

    Inverse Prob

    (2011)
  • Carlavan M, Blanc-Féraud L, Two constrained formulations for deblurring Poisson noisy images. In: Proc. IEEE...
  • T. Teuber et al.

    Minimization and parameter estimation for seminorm regularization models with I-divergence constraints

    Inverse Prob

    (2013)
  • H.W. Engl et al.

    Regularization of inverse problems

    (1996)
  • Cited by (0)

    This research is supported by the project Learning meets time: A new computational approach for learning in dynamic systems, contract RBFR12M3AC, funded by the Italian Ministry of Education, University and Research Futuro in Ricerca 2012 program and by the Italian Spinner2013 PhD project High-complexity inverse problems in biomedical applications and social systems.

    ☆☆

    The publisher would like to inform the readership that this article is a reprint of a previously published article. An error occurred which resulted in the publication of this article in a wrong issue. As a consequence, the publisher would like to make this reprint available for the reader’s convenience and for the continuity of the papers involved in the Special Issue. For citation purposes, please use the original publication details; Communications in Nonlinear Science and Numerical Simulation. Volume 20/3, March 2015, Pages 882-896. The publisher sincerely apologizes to the readership.

    View full text