A convergent blind deconvolution method for post-adaptive-optics astronomical imaging
aa r X i v : . [ a s t r o - ph . I M ] M a y A convergent blind deconvolution method forpost-adaptive-optics astronomical imaging
M Prato , A La Camera , S Bonettini and M Bertero Dipartimento di Scienze Fisiche, Informatiche e Matematiche, Universit`a di Modena eReggio Emilia, Via Campi 213/b, 41125 Modena, Italy Dipartimento di Informatica, Bioingegneria, Robotica e Ingegneria dei Sistemi, ViaDodecaneso 35, 16145 Genova, Italy Dipartimento di Matematica e Informatica, Universit`a di Ferrara, Via Saragat 1,44122 Ferrara, ItalyE-mail: [email protected] , [email protected] , [email protected] , [email protected] Abstract.
In this paper we propose a blind deconvolution method which applies todata perturbed by Poisson noise. The objective function is a generalized Kullback-Leibler(KL) divergence, depending on both the unknown object and unknown point spreadfunction (PSF), without the addition of regularization terms; constrained minimization,with suitable convex constraints on both unknowns, is considered. The problem is non-convex and we propose to solve it by means of an inexact alternating minimizationmethod, whose global convergence to stationary points of the objective function hasbeen recently proved in a general setting. The method is iterative and each iteration,also called outer iteration, consists of alternating an update of the object and the PSF bymeans of fixed numbers of iterations, also called inner iterations, of the scaled gradientprojection (SGP) method. Therefore the method is similar to other proposed methodsbased on the Richardson-Lucy (RL) algorithm, with SGP replacing RL. The use of SGPhas two advantages: first, it allows to prove global convergence of the blind method;secondly, it allows the introduction of different constraints on the object and the PSF.The specific constraint on the PSF, besides non-negativity and normalization, is anupper bound derived from the so-called Strehl ratio (SR), which is the ratio between thepeak value of an aberrated versus a perfect wavefront. Therefore a typical application,but not the unique one, is to the imaging of modern telescopes equipped with adaptiveoptics systems for partial correction of the aberrations due to atmospheric turbulence.In the paper we describe in detail the algorithm and we recall the results leading to itsconvergence. Moreover we illustrate its effectiveness by means of numerical experimentswhose results indicate that the method, pushed to convergence, is very promising in thereconstruction of non-dense stellar clusters. The case of more complex astronomicaltargets is also considered, but in this case regularization by early stopping of theouter iterations is required. However the proposed method, based on SGP, allowsthe generalization to the case of differentiable regularization terms added to the KLdivergence, even if this generalization is outside the scope of this paper. lind deconvolution in post-adaptive-optics astronomical imaging
1. Introduction
Blind deconvolution is the problem of image deblurring when the blur is unknown and,in general, is investigated by assuming a space-invariant model; in such a case the naiveproblem formulation is to solve the equation g = h ∗ f , where g is the detected image and ( f , h ) are respectively the unknown object and theunknown point spread function (PSF), while ∗ denotes convolution. It is obvious thatthis problem is extremely undetermined and that there is an infinite set of pairs solvingthe equation. Among them also the trivial solution ( f = g , h = δ ), where δ denotes theusual delta function. Therefore the problem must be reformulated by introducing as faras possible all available a priori information on both the object and the PSF.Blind deconvolution is the subject of a wide literature and we do not try to givea thorough account for that. Indeed the different approaches concern specific classesof images and PSFs. For instance approaches applicable to natural images may not besuitable in microscopy; approaches developed for motion blur are not applicable to otherclasses of blur, and so on. As concerns natural images we only mention a recent paper[32] which contains a critical analysis as well as several relevant references.In this paper we focus on astronomical imaging by assuming that an adaptiveoptics (AO) system is used to compensate for atmospheric blur and that a parametercharacteristic of this correction, the so-called Strehl ratio (SR), is approximately known.We recall that SR is the ratio of peak diffraction intensity of an aberrated versus perfectwaveform. In the case of AO images this parameter can be estimated by the astronomersduring the observation and provided with an error of few percent (about 4-5 %). Sincethis information provides an upper bound on the maximum value of the PSF, it can beused to exclude the trivial solution mentioned above and corresponding to the pair ( g , δ ).The approach we propose applies to astronomical imaging if noise is dominated byphoton counting and therefore the data are realizations of Poisson processes, even ifapproaches based on regularized least-square methods are also available (see for instance[29]). In the Poisson case several iterative methods have been already investigated, whichconsist of alternating updates of the object and PSF by means of Richardson-Lucy (RL)iterations [28, 23, 30, 42, 20, 21], or accelerated RL iterations [10].In [28] one iteration of the algorithm consists of updating both the object and the PSFby means of one RL iteration. This algorithm was investigated, in a different context, byLee and Seung [31] but their convergence proof is incomplete, since only the monotonicdecrease of the objective function is shown while, for a general descent method to beconvergent, strongest Armijo-like decreasing conditions have to be verified [35]. Theother approaches could be classified as methods of inexact alternating minimization sincethey use alternately a number of RL iterations on both the object and the PSF (remark, lind deconvolution in post-adaptive-optics astronomical imaging f and h : in the case of the object we only consider non-negativity while in the case of the PSF we consider both non-negativity and an upperbound provided by the knowledge of the SR, as well as the normalization condition whichmust be satisfied by the PSF. We point out that the relevance of the use of the SRconstraint for blind deconvolution was first pointed out by Desider`a & Carbillet [21] andthis paper intends to use it in a proper mathematical context.Thanks to [14] the convergence is assured if we use a fixed number of SGP iterationsfor updating the object and the PSF; the number of iterations may be different in the twocases (for the denomination asymmetric iterative blind deconvolution , see [10]). Since theproblem is non-convex, the limit of the iteration may depend on the choice of the initialstep and possibly on the numbers of internal iterations. The convergence result doesnot assure that the limit is a sensible solution of the problem, since we do not introduceregularization in our approach. A comment on this point is required.In the case of deconvolution it is well known that the minimizers of the discrepancyfunction for Poisson data, the generalized Kullback-Leibler (KL) divergence, are sparseobjects, i.e. they consist of bright spots over a black background; it is the so-called night-sky [1] or checker-board [34] effect. As a result these minimizers can be sensible solutionsin the case, for instance, of the deconvolution of images of not dense star clusters by a lind deconvolution in post-adaptive-optics astronomical imaging
2. Problem setting
Following [38], we assume that the observed image g can be modeled as the sum of twoterms g = g pe + r . The first, g pe , is the number of photo-electrons due to object andbackground emission and is a realization of a Poisson random variable with expected value ˜ g = ˜ h ∗ ˜ f + b , where ˜ f is the original object, ˜ h is the PSF of the acquisition system and lind deconvolution in post-adaptive-optics astronomical imaging b is the background term, while r represents the read-out noise (RON). Here and in thefollowing we denote by bold letters N × N arrays whose pixels are indexed by a multi-index i = ( i , i ), i ∈ S . For simplicity we assume that the background is constant andknown. As concerns the RON, it is a realization of a Gaussian additive random variablewith a known variance σ . According to Snyder et al. [39], it can be approximated by aPoisson process with mean and variance being the same as σ if the constant term σ isadded to g . If we add σ also to the background and if, with an abuse of notation, wedenote again as g and b the modified image and background, then we can conclude that g ∼ Poiss( ˜ h ∗ ˜ f + b ) . As concerns the PSF, we assume that it is normalized to unit volume X i ∈ S ˜ h i = 1 (1)and that its maximum value, denoted by s , is knownmax i ∈ S ˜ h i = s. (2)The upper bound s can be obtained by computing the diffraction-limited PSF of theconsidered telescope and multiplying its peak value by the SR value provided by theastronomers, as discussed in the Introduction.The blind deconvolution problem consists in finding an approximation of both ˜ f and ˜ h , given g , b and s . To this purpose we consider a maximum-likelihood approach to theproblem of image deconvolution. Since the maximization of the likelihood, which dependson the unknown object and PSF, is equivalent to the minimization of a generalized KLdivergence, we propose to estimate these approximations by minimizing this function (seethe comments in the Introduction concerning regularization) while taking into account allthe available information, i.e. the non-negativity of both the PSF and the original objectand the constraints (1)-(2). The resulting optimization problem is the followingmin KL ( g , h ∗ f + b ) (3)s . t . f ≥ ≤ h ≤ s , X i ∈ S h i = 1where KL denotes the KL divergence of h ∗ f + b from g KL ( g , h ∗ f + b ) = X i ∈ S (cid:26) g i log g i ( h ∗ f ) i + b i + ( h ∗ f ) i + b i − g i (cid:27) . (4)Problem (3) is convex if restricted to f or h only, but is in general non-convex withrespect to the pair ( f , h ), thus leading to the possible presence of several local minima. lind deconvolution in post-adaptive-optics astronomical imaging ∇ f KL ( g , h ∗ f + b ) = − H T gh ∗ f + b ∇ h KL ( g , h ∗ f + b ) = F T − F T gh ∗ f + b ∇ KL ( g , h ∗ f + b ) = H T diag (cid:16) g ( h ∗ f + b ) (cid:17) H H T diag (cid:16) g ( h ∗ f + b ) (cid:17) F − K ( h , f ) F T diag (cid:16) g ( h ∗ f + b ) (cid:17) H − K ( h , f ) T F T diag (cid:16) g ( h ∗ f + b ) (cid:17) F where H and F are the Block Circulant with Circulant Blocks (BCCB) matrices associatedto the convolution, i.e. h ∗ f = H f = F h , while K ( h , f ) is the Block Hankel with HankelBlocks (BHHB) matrix whose last row is the vector gh ∗ f + b (see [27, Chapter 4] for a surveyon structured matrices). Here the ratios and the squares are computed element-wise, and is a column vector with all entries equal to 1. Even if the diagonal blocks of the Hessianare symmetric positive semi-definite, the whole matrix is difficult to analyze and compute.
3. Alternating Minimization
Despite the complexity of the Hessian, the constraints have a simple, separable structure,which can be exploited by adopting an Alternating Minimization (AM) algorithm for thesolution of the non-convex problem (3)-(4). More precisely, the AM algorithms can beapplied to any problem of the formmin J ( x ) (5)s . t . x ∈ Ω × Ω × ... × Ω m ⊆ R n where, for all i = 1 , ..., m , Ω i is a closed and convex subset of R n i with n + ... + n m = n and any vector in the feasible set can be partitioned into vector components as x =( x , x , ..., x m ) x i ∈ Ω i . Clearly, the blind deconvolution problem (3) is a special caseof (5) with m = 2, x = f and x = h .The basic idea of AM is the cyclic minimization of the objective function withrespect to one variable, updating its value for the next optimization steps: in particular,AM is often referred to as the Nonlinear Gauss-Seidel (GS) method , where the iterate x ( k +1) = ( x ( k +1)1 , ..., x ( k +1) m ) is computed such that for i = 1 , ..., m the block of variables x ( k +1) i is a solution of the sub-problemmin J ( x ( k +1)1 , ..., x ( k +1) i − , y , x ( k ) i +1 , ..., x ( k ) m ) . (6)s . t . y ∈ Ω i This kind of approach has been widely studied in the literature [8, 9, 25, 26, 33, 41] andwe recall two important facts about it: lind deconvolution in post-adaptive-optics astronomical imaging • for m = 2 it has been proved in [26, Corollary 2] that the limit points of the sequence { x ( k ) } defined in (6) are stationary for problem (5) even in the non-convex case; • for m ≥ J : indeed, in [36], Powell devises a counterexample with m = 3 where all the limitpoints of the sequence generated by the nonlinear GS method are not stationary forthe problem (5). Some convergence results are proved for example in [8, 9, 26, 33]under suitable strict convexity assumptions.All the convergence results mentioned above, even in the case m = 2, are proved when theiterates are updated by an exact solution of the partial minimization problem (6), which isoften impractical or too costly to compute. Indeed, many practical AM algorithms, whichare also referred to as Block Coordinate Descent methods , are obtained by applying aniterative minimization method to approximately solve (6). In this case, the convergenceproperties of the alternating scheme also depend on the features of the inner solver. Adetailed analysis of the Block Coordinate Descent algorithms in the unconstrained caseis proposed in [25], where the authors devise some convergence conditions not necessarilyrelated to the convexity of the objective function.In this paper, we follow the approach in [14], where the partial minimization overeach variable (6) is performed inexactly by means of a fixed number of Scaled GradientProjection (SGP) steps [11].Choose a feasible starting point x (0) and a positive integer L ≥ k = 0 , , , ... $ For i = 1 , ..., m j Compute x ( k +1) i by applying n ( k +1) i ≤ L SGP iterations to (6) (7)A representation of the scheme (7) applied to the blind deconvolution problem isgiven in Algorithm 1: each main cycle consists of two successive deconvolution steps toupdate the current estimates of the object f ( k ) and PSF h ( k ) . For sake of completeness,we report the convergence result shown in Theorem 4.2 of [14] for our particular case. Theorem 3.1
Every limit point of the sequence ( f ( k ) , h ( k ) ) generated by Algorithm 1 isa stationary point for problem (3). As far as we know, convergence results stronger than that given in this Theorem forfirst-order methods applied to a general non-convex problem do not exist in the literature.The main difficulty in this kind of problems is to prove the existence of convergent sub-sequences. However, in our specific case, thanks to the Strehl constraint, the sequenceof the PSFs generated by our approach is bounded. Moreover, as concerns objectreconstruction one could introduce a constraint on the flux ( ℓ norm) of the reconstructed lind deconvolution in post-adaptive-optics astronomical imaging f , h ), the presence of multiple potential stationary points makes any limit pointdependent on both the initial guess and the chosen inner iteration numbers on the image( n f ) and the PSF ( n h ). Algorithm 1
Cyclic Scaled Gradient Projection (CSGP) methodChoose the starting point f (0) , h (0) and the inner iterations numbers n f , n h ≥ For k = 0 , , , ... do the following steps:Step 1. Compute f ( k +1) with n f SGP iterations applied tomin KL ( g , h ( k ) ∗ f + b ) (8)s . t . f ≥ f ( k ) Step 2.
Compute h ( k +1) with n h SGP iterations applied tomin KL ( g , h ∗ f ( k +1) + b ) (9)s . t . ≤ h ≤ s , X i ∈ S h i = 1starting from the point h ( k ) . End
We stress the fact that the main strength of Algorithm 1 for blind deconvolution withrespect to the more standard AM approach described for example in [16], is that it allowsan inexact solution of the inner sub-problems (8)–(9) while preserving the theoreticalconvergence properties. Since the proposed method is essentially based on SGP, we recallits main features in the following sub-section. lind deconvolution in post-adaptive-optics astronomical imaging The SGP algorithm is a first–order method which applies to any optimization problem ofthe form min x ∈ Ω J ( x ) (10)where J ( x ) is a continuously differentiable function and Ω is a convex set. Each SGPiteration is based on the feasible descent direction defined as d ( k ) = P Ω ,D − k ( x ( k ) − α k D k ∇ J ( x ( k ) )) − x ( k ) where α k is a scalar step-size parameter, D k is a diagonal matrix with positive diagonalentries and P Ω ,D − k ( · ) is the projection onto Ω associated to the norm induced by D − k , i.e. P Ω ,D − k ( x ) = arg min y ∈ Ω ( x − y ) T D − k ( x − y ) . (11)The new point is computed along the direction d ( k ) as follows x ( k +1) = x ( k ) + λ k d ( k ) where λ k is a step-length parameter to be chosen such that the monotone Armijo condition J ( x ( k ) + λ k d ( k ) ) ≤ J ( x ( k ) ) + βλ k ∇ J ( x ( k ) ) T d ( k ) (12)is satisfied for a fixed value of the parameter β ∈ (0 , λ k is computed by a standard backtrackingcondition as λ k = θ m , where θ ∈ (0 ,
1) and m is the smallest integer such that (12) issatisfied.The convergence of the SGP scheme, which is outlined in Algorithm 2, can be provedwhen the step-size α k and the diagonal entries of D k are bounded above and away fromzero, i.e. α k ∈ [ α min , α max ] with 0 < α min < α max and D k is chosen in the set D of diagonal matrices whose diagonal entries have values between L and L , for giventhresholds 0 < L < L .The SGP algorithm has been recently applied in several image restoration problems(see e.g. [3, 12, 43]). Under standard assumptions, it can be proved [11] that the SGPalgorithm is well defined and any limit point of the sequence { x ( k ) } is a stationary pointof (10); if, in addition, J ( x ) is convex, any limit point is a minimum point. It is worthstressing that the convergence result holds for any choice of the scaling matrix D k ∈ D and the step-length α k ∈ [ α min , α max ]: this freedom of choice can be exploited in orderto improve the convergence speed. Indeed, it is well known that gradient methods canbe significantly accelerated by a clever choice of the step-length parameter α k : one of themore effective strategies are the Barzilai–Borwein (BB) rules proposed firstly in [2] forquadratic unconstrained programming and then developed and analyzed for more generalproblems (see [17, 19, 24, 44] and reference therein). The BB rules can be considered a very lind deconvolution in post-adaptive-optics astronomical imaging Algorithm 2
Scaled gradient projection (SGP) methodChoose the starting point x (0) ≥ β, θ ∈ (0 , < α min < α max . For k = 0 , , , ... do the following steps:Step 1. Choose the parameter α k ∈ [ α min , α max ] and the scaling matrix D k ∈ D ; Step 2.
Compute the descent direction: d ( k ) = P Ω ,D − k ( x ( k ) − α k D k ∇ J ( x ( k ) )) − x ( k ) ; Step 3.
Backtracking loop: compute the smallest integer m such that (12) issatisfied with λ k = θ m ; Step 4.
Set x ( k +1) = x ( k ) + λ k d ( k ) . End cheap way to capture the second order information enforcing a quasi-Newton property.In our algorithm we adopt the scaled versions of the BB rules proposed in [11], which aregiven by α ( BB k = s ( k − T D − k D − k s ( k − s ( k − T D − k z ( k − , α ( BB k = s ( k − T D k z ( k − z ( k − T D k D k z ( k − , (13)where s ( k − = x ( k ) − x ( k − and z ( k − = ∇ J ( x ( k ) ) − ∇ J ( x ( k − ). Based on the previousformulas, we define the values α (1) k , α (2) k ∈ [ α min , α max ] in the following way if s ( k − T D − k z ( k − ≤ then α (1) k = min { · α k − , α max } ; else α (1) k = min n α max , max n α min , α ( BB k oo ; endifif s ( k − T D k z ( k − ≤ then α (2) k = min { · α k − , α max } ; else α (2) k = min n α max , max n α min , α ( BB k oo ; endif From the numerical experience, the best performances are obtained by an alternation ofthe two BB formulas: thus, following [24], we choose the following criterion for computing α k if k ≤ then α k = min j =max { ,k +1 − M α } ,...,k α (2) j ; (14) else if α (2) k /α (1) k ≤ τ k then Set α k as in (14) lind deconvolution in post-adaptive-optics astronomical imaging τ k +1 = 0 . · τ k ; else α k = α (1) k ; τ k +1 = 1 . · τ k ; endif where M α is a prefixed positive integer and τ ∈ (0 , D k = diag (cid:0) min( L , max( L , x ( k ) ) (cid:1) as suggested also in [11]. Obviously, in the inner iterations of (8), x is replaced by f while, in the inner iterations of (9), x is replaced by h . As concerns the choice of thebounds ( L , L ), at the beginning of each inner subproblem we perform one step of theRL method and tune the parameters according to the min/max positive values y min / y max of the resulting image according to the rule if y max /y min < then L = y min / L = y max · else L = y min ; L = y max ; endif Since the minimization steps (8) and (9) involve different constraints, corresponding totwo convex sets Ω and Ω , respectively, we have to account for two different algorithms tocompute the projections P Ω ,D − k and P Ω ,D − k . In the alternating procedure of Algorithm1, when SGP is applied to problem (8), the projection consists of a simple componentthresholding, obtained by setting all the negative elements of the vector to be projectedequal to zero. For the updating of the PSF, instead, we have to project on the constraintsset of the problem (9), consisting of a single linear equality constraint, in addition tosimple bounds (box constraints) on the variables. The resulting constrained optimizationproblem to be addressed is thereforemin 12 z T D − k z − z T y (15)s . t . ≤ z ≤ s , X i ∈ S z i = 1 lind deconvolution in post-adaptive-optics astronomical imaging y = D − k ( x ( k ) − α k D k ∇ J ( x ( k ) )). By introducing the Lagrangian penalty function,one can see that the orthogonal projection (15) can be re-conducted to a root-findingproblem of the piecewise linear monotonically non-decreasing function t ( ξ ) = X i ∈ S z i ( ξ ) − , where ξ is the Lagrangian multiplier of the equality constraint, z i ( ξ ) = mid(0 , ( D k ) ii ( y i + ξ ) , s )and mid( a , a , a §
4. Numerical experiments
In this section we investigate the effectiveness of the proposed blind method by means ofseveral numerical experiments. Since the blind problem formulated in (3) is non-convex,several local minima may exist. Moreover, we know that any limit point (or the limit)of the proposed iteration is a stationary point of the problem. The limit depends, ingeneral, on the numbers of inner iterations but also on the initialization of the outeriteration; therefore it is important to initialize the procedure with a sensible initial guessand we first discuss this point.As concerns the object we can use the standard initialization of the RL algorithm,namely a constant object with a flux coinciding with the flux of the image after backgroundsubtraction. The choice of the initial PSF is more important because, in the first step ofthe procedure, the image is deconvolved with this PSF.To this purpose we point out an important property of the PSF of a telescope: it is aband-limited function and, if the telescope consists of a circular mirror, the band, i.e. thesupport of its Fourier transform, is a disc with a radius proportional to the ratio betweenthe diameter D of the telescope and the observation wavelength λ . It is not easy to insertthis property as a constraint on the PSF because the projection on the resulting set ofconstraints (including SR and normalization) is not easily computable. For this reason wedo not consider this constraint in this paper. However we can try to force the estimatedPSF to have this property using an initial PSF which is band-limited and satisfies theother constraints.The ideal PSF of the telescope is not suitable as initial guess because it does notsatisfy the SR constraint. However one can consider, as suggested for instance in [10], lind deconvolution in post-adaptive-optics astronomical imaging ≥ s depends on both SR and the ratio D/λ ); for lowervalues of SR one can take the autocorrelation of the autocorrelation and so on, until theSR constraint is satisfied. This is the choice considered in our numerical experiments and,quite surprisingly, it seems that the algorithm, in spite of its high nonlinearity, preservesthe band-limiting property satisfied by the initial guess.All the numerical experiments have been performed with a set of routinesimplemented by ourselves in Interactive Data Language (IDL). The codes of the algorithmspresented and discussed in this paper are available under request.
As mentioned in the Introduction the use of SR as a constraint on the PSF is first proposedin [21]. Therefore some of our numerical experiments coincide with some of the testsperformed in that paper. In particular we use three of the AO-corrected PSFs (withSR equal to 0.67, 0.40 and 0.17, respectively), used by these authors and obtained bymeans of the Software Package CAOS [15]; the parameters corresponding to these PSFsare given in [21]. We only specify that they correspond to a telescope with an effectivediameter of 8.22 m and an observation wavelength of λ = 1 . µ m (H-band). For eachPSF, the images are generated by assuming, as in [21] a time exposure of 1200 s, with atotal transmission of 0.3. Moreover, a background of 13.5 mag arcsec − , corresponding toobservations in H-band, is added to the blurred images (for the convenience of the readerwe remark that it corresponds to 3 . × counts per pixel). The results are perturbedwith Poisson noise and additive Gaussian noise with σ = 10 e − / px. According to theapproach proposed in [39] and discussed in Sect. 2, RON compensation is obtained inthe deconvolution algorithms by adding the constant σ = 100 to the images and thebackground.In our first experiment we also consider an example which is not related to AOimaging but is a simulation of HST image before COSTAR correction, since this imageis frequently used in the testing of deconvolution methods. Obviously in such a case theideal PSF must be computed, by taking into account the diameter of the Hubble telescope,about 2.4 m, and the assumed observation wavelength of about 0.55 µ m, and compared tothe aberrated PSF in order to estimate the corresponding SR. Objects, PSFs and blurredimages used in all our experiments are sized 256 ×
256 pixels. lind deconvolution in post-adaptive-optics astronomical imaging We first report results on the following examples: • the binary system considered by Desider`a & Carbillet [21], in which the twocomponents have the same magnitude 12 in H-band (corresponding approximatelyto 6 . × counts) with an angular separation of 285 mas (19 pixels), i.e. ∼ ∼
40 mas); • a model of an open star cluster based on an image of the Pleiades, consisting of 9stars with magnitudes ranging from about 13 (i.e. about 2 . × counts) to 16(about 1 . × counts) in H-band and described in [7]; • a simulation of a star cluster, consisting of 470 light sources, as observed by theHubble Space Telescope (HST) before COSTAR correction. For this case only, wedo not use an AO-corrected PSF but the aberrated HST PSF, which corresponds toSR=0.09. These data can be obtained via anonymous ftp fromftp://ftp.stsci.edu/software/stsdas/testdata/restore/sims/star cluster/.In Fig. 1 we show the images of the binary and of the star cluster in the case of a PSFwith SR=0.67 as well as the HST image of a simulated star cluster. For all these exampleswe use 50 inner iterations on the object and 1 inner iteration for the PSF. This choice canbe justified by the features of our blind problem discussed in the Introduction, since weneed a sufficiently large number of SGP inner iterations for obtaining a nearly point-wiseobject. Moreover, with a few experiments on the binary, we verify that this choice is agood compromise which provides a sufficiently fast convergence for all cases. In a firstinstance we perform 300 outer iterations.As concerns the measure of the quality of the reconstructions, for the PSFs we use therelative r.m.s. error between the reconstructed PSF h and that used for image generation ˜ h , i.e. RM SE = || h − ˜ h || || ˜ h || , (16)where || . || denotes the ℓ norm. The same parameter is used for measuring the qualityof the reconstruction of the HST star cluster. In the case of the binary and of the openstar cluster we use a magnitude average relative error (MARE) defined as follows M ARE = 1 q q X i =1 | m i − ˜ m i | ˜ m i , (17)where q is the number of stars and m i , ˜ m i are respectively the reconstructed and the truemagnitudes.The results are shown in Table 1 and are consistent with the results reported in [21]but obtained with a sound mathematical approach, allowing investigation of the limit lind deconvolution in post-adaptive-optics astronomical imaging Figure 1.
Images of the binary and of the star cluster (upper panels) in the case of aPSF with SR=0.67. In the lower panel the image of the HST star cluster. for large number of iterations and generalization to regularized problems. The valuesof MARE estimated with our blind approach are certainly higher than those achievableif one deconvolves the data with the exact PSF ( inverse crime ) and given in the thirdcolumn, but they are still quite small. Moreover the reconstruction of the PSFs is verysatisfactory: the RMSEs of the initial PSFs are of the order of 30-50 %, while those ofthe reconstructed PSFs are of the order of few percents. A comparison between the true,initial and reconstructed PSFs is shown in Fig. 2. We must add that the reconstructionerror is still decreasing after 300 iterations and therefore the minimum of the objectivefunction is not yet reached.For investigating the limit of the algorithm we consider three other examples ofbinaries: a binary with magnitudes 12-12 and angular distance of 10 pixels and twobinaries with magnitudes 12-16 and angular distances of 19 and 10 pixels respectively.For these four examples of binaries we generate images using the PSF with the highestSR, namely 0.67, and we compute 8000 outer iterations of the blind algorithm, usingthe fixed pair ( n f , n h ) = (50 , lind deconvolution in post-adaptive-optics astronomical imaging Figure 2.
First column: the PSFs used for image generation; second column: thePSFs used for initializing the blind algorithm (see the text for their computation); thirdcolumn: the PSFs reconstructed by the blind algorithm. First row: AO-corrected PSFwith SR=0.67; second row: AO-corrected PSF with SR=0.40; third row: AO-correctedPSF with SR=0.17; fourth row: HST PSF before COSTAR correction. lind deconvolution in post-adaptive-optics astronomical imaging Table 1.
Reconstruction errors for point-wise objects. In the first column we specify theobject and in the second the value of the SR used for image generation; in the third andfourth the values of MARE (RMSE in the case of HST cluster) when SGP is used forimage deconvolution with the exact and initial PSF respectively. In the fifth column thevalues of MARE (RMSE in the case of HST) obtained with 300 outer iterations by ourblind approach, using the fixed pair ( n f , n h ) = (50 , Image SR
M ARE M ARE M ARE RM SE RM SE Binary 0.67 1 . × − . × − . × −
32 % 1 . . × − . × − . × −
54 % 2 . . × − . × − . × −
55 % 3 . . × − . × − . × −
32 % 1 . . × − . × − . × −
54 % 1 . . × − . × − . × −
55 % 4 . . . . × − .It should be interesting to find a way for establishing if the minima we find are theglobal ones or not, but, as it is known, global minimization is a very difficult problem.As a test, even if it does not provide a proof that the minima are the global ones, wecompare the minimum values of the objective function, i.e. the KL divergence, with itsvalues corresponding to the ground truths, i.e. the values obtained by substituting in Eq.(4) the objects and PSFs used for image generation. We find that the minimum values areof the order of 1 . × while the values corresponding to the ground truth are greater byabout a factor of 3. We can only say that, if these values were smaller than our minimum lind deconvolution in post-adaptive-optics astronomical imaging Table 2.
Reconstruction errors, provided by increasing number of iterations, in the caseof four different binaries (the parameters are indicated in the first column, as explainedin the text) whose images are generated using the PSF with SR=0.67. As usual MAREis a measure of the errors on the magnitudes of the two stars while RMSE is a measureof the error on the reconstructed PSF.
300 it 4000 it 8000 it12-12 MARE 1 . × − . × − . × −
19 pixels RMSE 1 . .
17 % 0 .
15 %12-16 MARE 5 . × − . × − . × −
19 pixels RMSE 0 .
99 % 0 .
20 % 0 .
27 %12-12 MARE 1 . × − . × − . × −
10 pixels RMSE 1 . .
18 % 0 .
17 %12-16 MARE 1 . × − . × − . × −
10 pixels RMSE 1 . . . m = 8. For this reason we performan analysis of this problem in the case of the binary, using the PSF with the highest SR,namely SR=0.67. We first consider the inverse crime reconstruction of the binary with magnitudes 12-12and angular distance 19 pixels. We deconvolve its image generated by means of the PSFwith SR=0.67 using the same PSF (inverse crime). The algorithm is the standard SGPwith non-negativity constraint. Also in this case the reconstruction is not free of artifactsbut they are randomly distributed and their values are quite small: the brightest artifacthas a magnitude m = 24, hence with a difference ∆ m = 12 with respect to the stars ofthe binary.Next we apply the blind algorithm, using as a constraint the exact value of SR andwe analyze the results provided by the 8000 outer iterations already considered in theprevious section. In the first two rows of Fig. 3 we show the reconstructed PSF and the lind deconvolution in post-adaptive-optics astronomical imaging RM SE = 1 . RM SE = 0 .
17 %
RM SE = 0 .
15 %
M ARE = 1 . × − M ARE = 1 . × − M ARE = 1 . × − RM SE = 2 . RM SE = 2 . RM SE = 2 . M ARE = 8 . × − M ARE = 4 . × − M ARE = 4 . × − Figure 3.
Binary with magnitudes 12-12 and angular separation of 19 pixels. First andsecond row: the reconstructed PSF and object after 300 (left), 4000 (middle), and 8000(right) outer iterations with the exact value of SR (SR=0.67) as a constraint. Third andfourth rows: the reconstructed PSF and object after 300 (left), 500 (middle), and 2000(right) outer iterations with the underestimated value of SR (SR=0.64) as a constraint.The symbol lind deconvolution in post-adaptive-optics astronomical imaging . m = 10 and after 8000 iterations ∆ m = 11. We find similar results in thecase of the open star cluster and therefore we can conclude that a very large number ofiterations may be required for obtaining very good results, at least if the exact value ofSR is known.However, as briefly discussed in the Introduction, it is not possible to know exactlythe value of SR. According to astronomers the expected error is of the order of 4 %.Therefore, we consider a variation of the constraint of this order of magnitude for imagesgenerated in the case of SR=0.67; more precisely we consider two values, SR=0 . .
64. We apply the blind algorithm using as a constraint the corresponding valuesof s . In the first case the reconstructions are definitely worse, the number of artifactsconsiderably increases as well as the error on the PSF. For instance, in the case 12-12 theRMSE is of the order of 5 % and does not decrease with increasing number of iterations(remember that, in the case of exact value, the error after 8000 iterations is of the order of0.15 %). On the other hand, if we underestimate the SR, i.e. we take as a constraint thevalue of s corresponding to SR=0.64, then the results are satisfactory. The reconstructionerrors for the four binaries already considered in the previous subsection are reported inTable 3. By comparing with the results reported in Table 2 and referring to the exactconstraint, we can conclude that the reconstruction errors are not significantly greaterthan those obtained in the exact case and that the convergence is faster.In the case of the binary with magnitudes 12-12 and angular distance 19 pixels weshow the reconstructions of the PSF and of the binary after 300, 500 and 2000 iterationsrespectively in the third and fourth row of Fig. 3. Quite surprisingly, the artifacts in thereconstruction of the binary completely disappear after 2000 iterations and the error onthe reconstructed PSF is of the order of 2 %. We can add that the same result is obtainedin the case of the open star cluster.In order to further investigate the effect of a wrong value of s , in the case of the lind deconvolution in post-adaptive-optics astronomical imaging Table 3.
Reconstruction errors in the case of the four different binaries of Table 2(described in the text), considering an underestimated constraint of the blind algorithm(SR=0.64). As usual MARE is a measure of the errors on the magnitudes of the twostars while RMSE is a measure of the error on the reconstructed PSF.
300 it 500 it 2000 it12-12 MARE 8 . × − . × − . × −
19 pixels RMSE 2 . . . . × − . × − . × −
19 pixels RMSE 2 . . . . × − . × − . × −
10 pixels RMSE 2 . . . . × − . × − —10 pixels RMSE 2 . . . As additional examples of astronomical targets, we consider three HST images: the Crabnebula NGC 1952, the galaxy NGC 6946 and the planetary nebula NGC 7027. In allcases we assume an integrated magnitude equal to 10 and, for each one, we obtain threeblurred images by convolving with the three PSFs of SR=0.67, 0.40 and 0.17. Again abackground in H-band is added to all images and the results are perturbed with Poissonand additive Gaussian noise. In the first column of Fig. 4 we show the three objects inreversed scale of gray levels while in the second column we show their blurred images inthe case SR=0.67.As concerns the initialization of the blind algorithm we use the same PSFs alreadyused in the case of stellar objects, namely obtained by suitable autocorrelations of theideal PSF of the telescope. However in the case of complex and diffuse objects, as alreadyremarked by other authors (see, for instance, [10]), a difficult and crucial point is thechoice of the number of inner iterations. We do not have a rule which can be successfullyapplied to all cases as for the stellar objects, i.e. ( n f , n h ) = (50 , lind deconvolution in post-adaptive-optics astronomical imaging Table 4.
Reconstruction errors for complex and diffuse objects. In third and fourthcolumns, the best errors achieved by SGP with the true and the initial PSFs, respectively.In the fifth column, the best error obtained using a maximum of 100 outer iterations (forthe choice of the inner iterations see the text). Finally, in the last two columns, the errorbetween the true PSF and the initial one, followed by the error between the true PSFand the one obtained in conjunction with the best reconstruction of the correspondingobject.
Image SR
RM SE obj
RM SE obj RM SE obj RM SE psf RM SE psf Crab 0.67 11 % 16 % 12 % 32 % 6 . . . . . . . . . . . n f , n h ) = (13 ,
22) forSR=0.67 and ( n f , n h ) = (13 ,
27) for SR=0.40, if we search for minimum RMSE on theobject using 100 outer iterations. Moreover, even if the algorithm is convergent, the limitis not, in general, a sensible solution: a suitable stopping of the outer iterations is required.In other words the algorithm is semi-convergent [34, 4] as concerns the outer iterations,i.e. the RMSE on the object first decreases, reaches a minimum and then goes away. Wedo not have a proof of this feature which derives from the numerical experiments. It isobvious that such a situation is not satisfactory and we briefly discuss this point in thenext section (see also the Introduction).In Table 4 we report the best results we have obtained while in the third column ofFig. 4 we show the reconstructions of the objects provided by the blind algorithm. Westress the case of the planetary nebula: it seems that in this case the algorithm is unableto improve the reconstruction with respect to that provided by the initial guess. We alsoremark that the error on the reconstructed PSF depends on the object: for instance, forthe galaxy it is larger than that for the Crab.We conclude by reporting an experiment intended to check the effect of anunderestimated or overestimated SR in the reconstruction of a diffuse object. We considerthe case of the Crab nebula and PSF with SR=0.67, and we apply our algorithm withSR=0.6 and SR=0.8. In the first case, the minimum RMSE on the object and the PSF areequal to 13% and 8.2%, while in the other the errors on object and PSF are 12% and 6.7%, lind deconvolution in post-adaptive-optics astronomical imaging Figure 4.
First column: the objects used for image generation; second column: theblurred images in the case of SR=0.67; third column the reconstructed objects obtainedwith our blind algorithm. First row: the Crab nebula NGC 1952, second row: the spiralgalaxy NGC 6946, third row: the planetary nebula NGC 7027. which are essentially the same values obtained with the correct SR. The small varianceobserved suggests that, in presence of complex objects, the availability of a correct SRdoes not represent a crucial point. In this case, an overestimate of the SR value seems tobe preferable.
5. Concluding remarks and perspectives
In this paper we propose a blind algorithm, based on SGP, for the reconstruction ofastronomical images acquired by a telescope equipped with an AO system. The algorithmcan be classified as an inexact alternating minimization of the KL divergence depending on lind deconvolution in post-adaptive-optics astronomical imaging
StarFinder [22], developed for accurate photometric and astrometric analysisof star clusters. These codes contain a method for extracting the PSF from the image ofthe star field; this PSF can be compared with and/or replaced by that provided by ourblind approach; also in this case the analysis of simulated star fields, in particular crowdystar fields, can help to understand when the blind approach is required.As already remarked in the Introduction, the situation is not so clear in the caseof more complex astronomical targets. Our preliminary simulations indicate that theproposed blind method has the semi-convergent property as concerns the outer iterations(the numbers of the inner iterations are fixed by the user and, in any case, they shouldnot be too large). The difficulties found in optimizing the numbers of inner iterationsmay reside in the fact that in this paper we do not introduce regularization in theobjective function. Due to the flexibility of SGP the method can be easily generalized todifferentiable regularizers (one only needs to compute the gradient of the penalty and itspositive part) and this generalization will be the subject of future work. We stress againthat the crucial point is to identify regularizers which are suitable for specific classes ofastronomical targets as well as regularizers which are suitable for AO corrected PSFs.The codes of the algorithms presented and used in this paper are available underrequest.
Acknowledgments
This work has been partially supported by MIUR (Italian Ministry for Universityand Research), PRIN2008 “Optimization Methods and Software for Inverse Problems”,grant 2008T5KA4L, and FIRB - Futuro in Ricerca 2012 “Learning meets time: a newcomputational approach for learning in dynamic systems”, contract RBFR12M3AC 002,and by INAF (National Institute for Astrophysics) under the contract TECNO-INAF2010 “Exploiting the adaptive power: a dedicated free software to optimize and maximize lind deconvolution in post-adaptive-optics astronomical imaging
References [1] Barrett H H and Meyers K J 2003
Foundations of Image Science (New York: Wiley and Sons)[2] Barzilai J and Borwein J M 1988 Two point step size gradient methods
IMA J. Numer. Anal. Inverse Problems Introduction to Inverse Problems in Imaging (Bristol: IoPPublishing)[5] Bertero M, Bindi D, Boccacci P, Cattaneo M, EVA C and Lanza 1998 A novel blind deconvolutionmethod with an application to seismology
Inverse Problems Inverse Problems Inverse Problems Nonlinear Programming (Belmont: Athena Scientific)[9] Bertsekas D P and Tsitsiklis J 1988
Parallel and Distributed Computation: Numerical Methods (Belmont: Prentice-Hall)[10] Biggs D S C and Andrews M 1998 Asymmetric iterative blind deconvolution of multi-frame images
Proc. SPIE
Inverse Problems Inverse Problems Inverse Problems IMA J. Numer. Anal. Mon. Not. R. Astron. Soc.
Linear Algebra Appl.
Math.Program.
Math. Program.
IMA J. Numer. Anal. Astron. Astrophys. lind deconvolution in post-adaptive-optics astronomical imaging [21] Desider`a G and Carbillet M 2009 Strehl-constrained iterative blind deconvolution for post-adaptive-optics data Astron. Astrophys.
Proc.SPIE
J. Opt. Soc. Am.
A-12
J. Ind. Manag. Optim. Optim. Method Softw. Oper. Res. Lett. Deblurring Images. Matrices, Spectra, and Filtering (Philadelphia: SIAM)[28] Holmes T J 1992 Blind deconvolution of quantum-limited incoherent imagery: maximum-likelihoodapproach
J. Opt. Soc. Am.
A-9
Astrophys. J.
Proc. SPIE
Advances in NeuralInformation Processing 13 (Proc. NIPS ∗ (MIT Press) 556–62[32] Levin A Weiss Y Durand F and Freeman W T 2009 Understanding and evaluating blinddeconvolution algorithms IEEE Conf. Computer Vision and Pattern Recognition
J. Optimiz. Theory App. Mathematical Methods in Image Reconstruction (Philadelphia:SIAM)[35] Nocedal J and Wright S J 2006
Numerical Optimization: Second Edition (New York: Springer)[36] Powell M J D 1973 On search directions for minimization algorithms
Math. Program. Astron. Astrophys.
A133[38] Snyder D L, Hammoud A M and White R L 1993 Image recovery from data acquired with acharge-coupled-device camera
J. Opt. Soc. Am.
A-10
J. Opt. Soc. Am.
A-12
Inverse Problems J. Optimiz. TheoryApp. Astron. Astrophys.
Inverse Problems Comput. Optim. Appl.35