ACQUIRE: an inexact iteratively reweighted norm approach for TV-based Poisson image restoration
Daniela di Serafino Germana Landi Marco Viola

	





[t1]This work
	was partially supported by Gruppo Nazionale per il Calcolo Scientifico - Istituto Nazionale di Alta Matematica (GNCS-INdAM).


1]Daniela di Serafinocor1daniela.diserafino@unicampania.it
2]Germana Landigermana.landi@unibo.it
3]Marco Violamarco.viola@uniroma1.it

[1]Dipartimento di Matematica e Fisica, Università degli Studi della Campania 
 “Luigi Vanvitelli”, viale A. Lincoln 5, 81100 Caserta, Italy
[2]Dipartimento di Matematica, Università degli Studi di Bologna, Piazza di Porta S. Donato 5, 40126 Bologna, Italy
[3]Dipartimento di Ingegneria Informatica, Automatica e Gestionale “Antonio Ruberti”, Sapienza Università di Roma, via Ariosto 25, 00185 Roma, Italy


We propose a method, called ACQUIRE, for the solution of constrained optimization problems modeling the restoration of images
corrupted by Poisson noise. The objective function is the sum of a generalized Kullback-Leibler divergence term and a TV regularizer,
subject to nonnegativity and possibly other constraints, such as flux conservation. ACQUIRE is a line-search method that considers
a smoothed version of TV, based on a Huber-like function, and computes the search directions by minimizing
quadratic approximations of the problem, built by exploiting some second-order information. A classical second-order Taylor
approximation is used for the Kullback-Leibler term and an iteratively reweighted norm approach for the smoothed TV term.
We prove that the sequence generated by the method has a subsequence converging to a minimizer of the smoothed problem
and any limit point is a minimizer. Furthermore, if the problem is strictly convex, the whole sequence is convergent. We note
that convergence is achieved without requiring the exact minimization of the quadratic subproblems; low accuracy in this
minimization can be used in practice, as shown by numerical results. Experiments on reference test problems show that our method
is competitive with well-established methods for TV-based Poisson image restoration, in terms of both computational
efficiency and image quality.


image restoration, Poisson noise, TV regularization, iteratively reweighted norm approach, quadratic approximation. 








§ INTRODUCTION


Restoring images corrupted by Poisson noise is required in many applications, such as fluorescence microscopy <cit.>,
X-ray computed tomography (CT) <cit.>, positron emission tomography (PET) <cit.>,
confocal microscopy <cit.> and astronomical imaging <cit.>.
Thus, this is a very active research area in image processing. We consider a discrete formulation of the problem, where the object
to be restored is represented by a vector
∈^n and the measured data are assumed to be a vector ∈_0^m, whose entries y_j
are samples from m independent Poisson random variables Y_j with probability

    P(Y_j = y_j) = e^-(A  + )_j(A + )_ j^y_j/y_j!,

where the matrix A = (a_ij) ∈^m × n models the observation mechanism of the imaging system
and ∈^m,  > 0, models the background radiation detected by the sensors.
Standard assumptions on A are

    a_ij≥ 0   i, j,     ∑_i=1^m a_ij = 1  j.

By applying a maximum-likelihood approach <cit.>,
we can estimate  by minimizing the Kullback-Leibler (KL) divergence of A + from :

    D_KL(A  + , ) = ∑_j=1^m ( y_j lny_j/(A +)_j + (A +)_j - y_j ),

where we set y_j ln (y_j / (A +)_j) = 0 if y_j = 0.
A regularization term is usually added to (<ref>) to deal with the inherent ill-conditioning of the estimation problem.
We focus on edge-preserving regularization by Total Variation (TV) <cit.>,
which has received considerable attention because of
its ability of preserving edges and smoothing flat areas of the images. We note that, although TV regularization
is known to suffer from undesirable staircase artifacts, it is still widely used in many medical and biological applications (see, e.g., <cit.>, <http://ranger.uta.edu/ huang/R_CSMRI.htm>).
Furthermore, by focusing on TV-regularized problems, we introduce and test a novel solution method that allows for extensions
to other models, such as high-order TV <cit.> and Total Generalized
Variation <cit.>, proposed to reduce the staircase effect.

Assuming, for simplicity, that  is obtained by stacking the columns of a 2D image X = (X_k,l) ∈^r × s, i.e.,
x_i = X_k,l with i = (l-1)r + k and n=rs, the following discrete version of the TV functional can be defined <cit.>:

    TV(X) = ∑_k=1^r∑_l=1^s√((X_k+1,l-X_k,l)^2 + (X_k,l+1-X_k,l)^2),

where X is supposed to satisfy some boundary conditions, e.g., periodic. This can be also
written as

    TV() = ∑_i=1^n D_i ,

where

    D_i =
     ( [ _(l-1)r + k + 1^T - _(l-1)r + k^T;           _lr+k^T - _(l-1)r + k^T ]),
       
    i = (l-1)r + k,

_q ∈^n is the qth standard basis vector, and · is the 2-norm


















Thus, we are interested in solving the following problem:

    [                    D_KL()  + λ  TV(),;                                    ∈, ]

where D_KL() is a shorthand for D_KL (A  + , ), λ > 0
is a regularization parameter, and ∈ corresponds to some physical constraints.
The nonnegativity of the image intensity naturally leads to the constraint ≥0̱.
When the matrix A comes from the discretization of a convolution operator and it is normalized as
in (<ref>), the constraint ∑_i=1^n x_i = ∑_j=1^m (y_i-b_i)
can be added, since the convolution performs a modification of the intensity distribution, while the total intensity remains
constants <cit.>.[We have implicitly assumed that  has been converted into
a real vector with entries ranging in the same interval as the entries of .]
In other words, common choices of  are

    = _1 := {∈^n : ≥0̱}

or

    = _2 :=  {∈^n : ≥0̱,  ^T  =  ^T ( - ) },

where  and  denote the vectors of all 1's of sizes n and m, respectively.

Various approaches have been proposed to solve problem (<ref>), mostly with = _1;
a key issue in all cases is to deal with the nondifferentiability of the TV functional.
Some representative methods are listed next.
A classical approach consists in approximating the TV functional with a smooth version of it
and using well-established techniques such as expectation-maximization methods
<cit.>,
gradient-projection methods with suitable scaling techniques aimed at accelerating convergence
<cit.>,
and alternating linearized minimization methods <cit.>.

The approximation of TV can be avoided, e.g., by using forward-backward splitting techniques;
this is the case of the proximal-gradient methods proposed in <cit.>
and the forward-backward EM method discussed in <cit.>.
On the other hand, the previous methods require, at each step, the solution of a Rudin-Osher-Fatemi (ROF) denoising
subproblem <cit.>, which can be computed only approximately, using, e.g., the algorithms
proposed in <cit.>.
Methods based on ADMM and SPLIT BREGMAN techniques, such as those presented in
<cit.>, do not exploit smooth TV
approximations too. They generally use more memory because of auxiliary variables of the same size
as  or , and require the solution of linear systems involving A^TA and, possibly,
the solution of ROF subproblems. Finally, a different approach to avoid the difficulties associated with
the nondifferentiability of the TV functional is based on the idea of reformulating (<ref>) as
a saddle-point problem and solving it by a primal-dual algorithm. In this context, an alternating extragradient
scheme has been presented in <cit.>, and a procedure exploiting the
Chambolle-Pock algorithm <cit.> has been described in <cit.>.












In this paper we take a different approach, aimed at exploiting some second-order information not considered
by the aforementioned methods. We consider a smoothed version of TV, based on a Huber-like function,
and propose a line-search method, called ACQUIRE (Algorithm based on Consecutive QUadratic and Iteratively
REweighted norm approximations), which minimizes a sequence of quadratic models obtained by
a second-order Taylor approximation of the KL divergence and an iteratively reweighted norm (IRN) approximation of the smoothed TV.
We prove the convergence of ACQUIRE with inexact solution of the inner quadratic problems. We show by numerical
experiments that exploiting some second-order information can lead to fast image restorations even with
low accuracy requirements on the solution of the inner problems, without affecting the quality of the reconstructed images.
In particular, ACQUIRE generally produces a strong reduction of the reconstruction error in the first iterations,
thus achieving a good tradeoff between accuracy and efficiency, and resulting competitive with
state-of-the-art methods.

The remainder of this paper is organized as follows. In Section <ref> we recall some preliminary
concepts that will be exploited later. In Section <ref> we describe our method and in Section <ref>
we prove that it is well posed and convergent. We provide implementation details and discuss the results obtained
by applying the proposed method to several test problems in Section <ref>. Some conclusions are reported in
Section <ref>.




§ PRELIMINARIES


We first provide some useful details about the KL divergence and introduce a smooth version of the TV functional.
Then we recall the concept of projected gradient and its basic properties, exploited later in this work.

Assumptions (<ref>) and > 0 ensure that, for any given ≥ 0, D_KL is a nonnegative,
convex, coercive, twice continuosly differentiable function in _+^n (see, e.g.,
<cit.>). Its gradient and Hessian
are given by

    ∇ D_KL()= A^T (  - /A  + )

and

    ∇^2 D_KL() = A^T U()^2 A,    U() = ( √()/A  + ),

where the square root and the ratios are intended componentwise, and () denotes the diagonal matrix with
diagonal entries equal to the entries of . It can be proved that ∇ D_KL is Lipschitz
continuous <cit.>; furthermore, it follows from (<ref>)
that ∇^2 D_KL is positive definite, i.e., D_KL is strictly convex, whenever > 0̱ and A has nullspace
(A) = {0̱}. In this case, if  is constrained to be in a bounded subset of the nonnegative orthant, e.g.,
the set _2 in (<ref>), the minimum eigenvalue of ∇^2 D_KL()  is bounded below
independently of , and D_KL is strongly convex.

From a practical point of view, it is interesting to note that A is usually the representation of a convolution operator, and hence
the computation of ∇ D_KL or of matrix-vector products involving ∇^2 D_KL
can be performed efficiently via fast algorithms for discrete Fourier, cosine or sine transforms.

The  TV functional is nonnegative, convex and continuous. Thus problem (<ref>) admits
a solution, which is unique if  > 0 and (A) = {0̱}. Since TV is not differentiable,
we use a regularized version of it, TV_μ. Taking into account the discussion
in <cit.> about smoothed versions of TV, we consider

    TV_μ() = ∑_i=1^nϕ_μ(D_i ),

where ϕ_μ is the Huber-like function

    ϕ_μ(z) = {[          |z|     |z| > μ,; 1/2(z^2/μ+μ)            , ].

and μ > 0 is small. It is easy to verify that TV_μ is Lipschitz continuously differentiable and its gradient reads as follows:

    ∇ TV_μ () = ∑_i=1^n∇ϕ_μ(D_i ),   ∇ϕ_μ( D_i ) = {[ D_i^T D_i / D_i       D_i  > μ,;    D_i^T D_i /μ               . ].

We also observe that TV_μ is not twice continuously differentiable, but has continuous Hessian for all
 such that D_i μ:

    ∇^2 TV_μ () = ∑_i=1^n∇^2 ϕ_μ(D_i ),


    ∇^2 ϕ_μ( D_i ) =  {[ D_i^T D_i/ D_i  - (D_i^T D_i ) (D_i^T D_i )^T/D_i ^3                                           D_i  > μ ,;                                          D_i^T D_i/μ                                           D_i  < μ . ] .



Now we recall basic notions about the projected gradient.
Let  be a nonempty, closed and convex set. For any continuously differentiable
function f: ⊆^n →, with  open set containing , the projected gradient
of f at ∈ is defined as the orthogonal projection of -∇ f onto the tangent cone
to  at , denoted by T_():

    ∇_ f () = min{‖ + ∇ f () ‖  ∈ T_() },

When  is the set _1 defined in (<ref>), the tangent cone takes the form

    T_() = {∈^n   :   v_i ≥ 0  x_i  = 0 }

and the computation of ∇_ f () is straightforward; when  is the set
_2 in (<ref>),

    T_ () = {∈^n   :  ^T=0  v_i ≥ 0  x_i  = 0 },

and ∇_ f () can be efficiently determined too, thanks to
the availability of low-cost algorithms for computing the projection in this case (see, e.g.,
<cit.>).

Since the projection onto  is nonexpansive, for all , ∈ it is

    ∇_ f () - ∇_ f () ≤ - ;

furthermore,

    - ∇ f() = ∇_ f() + P_N_() (- ∇ f()),

where P_N_() denotes the orthogonal projection operator onto the normal cone to  at ,

    N_ () = {∈^n   :  ^T ≤ 0 ∈() },

which is the polar cone of T_ () (see, e.g., <cit.>).


Finally, it is well known that any constrained stationary point ^* of f in  is characterized by
∇_ f (^*) = 0̱










and that ∇_ f is lower semicontinuous on  (see, e.g., <cit.>).




§ IRN-BASED INEXACT MINIMIZATION METHOD


We propose an iterative method for solving the problem

    [                       D_KL()  + λ  TV_μ (),;                                          ∈, ]

where  can be any nonempty, closed and convex subset of R_+^n, although our practical interest
is for the feasible sets in (<ref>)-(<ref>).
This method is based on two main steps: the inexact solution of a quadratic model of (<ref>)
and a line-search procedure.

Given an iterate ^(k)∈, we consider the following quadratic approximation of D_KL:

    [                                D_KL ()  ≈ D_KL^(k)()                = D_KL(^(k)) + (-^(k))^T ∇ D_KL(^(k));                                                      +   1/2 (-^(k))^T ( ∇^2 D_KL (^(k)) + γ I ) (-^(k)), ]

where I is the identity matrix and γ > 0. Note that γ I has been introduced to ensure that D_KL^(k) is strongly convex;
obviously, we can set γ = 0 if > 0̱, (A) = {0̱} and  is bounded.

In order to build a quadratic model of TV_μ, we use the IRN approach described in <cit.>,
i.e., we approximate TV_μ() as follows:

    TV_μ() ≈ TV_μ^(k)() = 1/2∑_i=1^n w_i^(k) D_i ^2 + 1/2 TV_μ(^(k)),

where








    w_i^(k) = {[   D_i ^(k)^-1 D_i ^(k) > μ,;          μ^-1             . ].

Trivially,

    TV_μ^(k) (^(k)) = TV_μ (^(k)),   ∇ TV_μ^(k) (^(k))  = ∇ TV_μ (^(k)) ;






furthermore,

    ∇^2 TV_μ^(k) (^(k)) = ∑_i=1^n w_i^(k) D_i^T D_i ,

and hence, for any  such that D_i ^(k)μ, the Hessian ∇^2 TV_μ^(k)(^(k)) can
be regarded as an approximation of ∇^2 TV_μ (^(k)), obtained by neglecting the higher order term
in the right-hand side of (<ref>), which generally increases the ill-conditioning of the Hessian matrix.




Thus, we can say that TV_μ^(k) contains some second-order information about TV_μ.
It is worth noting that the higher order term of the Hessian of a smoothed TV function is also neglected
in the lagged diffusivity method by Vogel and Oman <cit.>.

In the following, to simplify the notation we set

    F()             =     D_KL () + λ  TV_μ(), 
    
       F_k()     =     D_KL^(k) () + λ  TV_μ^(k) ().

At iteration k, our method computes a feasible approximation ^(k) to the solution ^(k)
of the quadratic problem

    [         F_k (),;              ∈, ]

and performs a line search along the direction

    ^(k) = ^(k) - ^(k),

until an Armijo condition is satisfied, to obtain an approximation ^(k+1) to the solution of problem (<ref>).
This procedure is sketched in Algorithm <ref> and is called ACQUIRE, which comes from “Algorithm based
on Consecutive QUadratic and Iteratively REweighted norm approximations”.







ACQUIRE is well posed (i.e., a steplength α_k satisfying the Armijo condition
can be found in a finite number of iterations) and is convergent; this is proved in Section <ref>.
Step <ref> does not require the exact solution of problem (<ref>), but only
the computation of an approximate solution such that condition (<ref>) at line <ref>
of the algorithm holds, with lim_k →∞ε_k = 0.

In Section <ref> we also show that the first condition in (<ref>) is satisfied if

    ∇_ F_k (^(k)) ≤θ^k  ∇_ F_k (^(0))  ,

and θ∈ (0,1). Therefore, the first condition in (<ref>) can be replaced by another one
which is simple to verify when the projected gradient can be easily computed, e.g.,
in the practical cases where  is one of the sets in (<ref>)-(<ref>).

The second condition in (<ref>) can be achieved by using any constrained minimization algorithm.
We note that, for the restoration problems considered in this work, gradient-projection methods, such as those
in <cit.>,
are suited to the solution of the inner problems (<ref>). Indeed, numerical experiments
have shown that very low accuracy is required in practice in the solution of the inner problems;
furthermore, the computational cost per iteration of gradient projection methods is modest
when low-cost algorithms for the projection onto the feasible set are available. More details on the
inner method used in our experiments are given in Section <ref>.




§ WELL-POSEDNESS AND CONVERGENCE


In order to prove that ACQUIRE is well posed, we need the following lemma <cit.>.


Let f : ^n → be continuously differentiable and let , ∈^n. If
there exists L > 0 such that

    ∇ f( + t ) - ∇ f () ≤ L t   t ∈ [0,1] ,

then

    f ( + ) ≤ f() + ∇ f ()^T  + L/2^2.



We also observe that, at step <ref> of Algorithm <ref>, we can find ^(k)^(k)
unless ^(k) is the solution ^(k) of problem (<ref>). However, in this case
^(k) is the solution of problem (<ref>), since the gradients, and hence the projected gradients,
of the objective functions of the two problems coincide at ^(k). Therefore, in the following
we can assume that ^(k)^(k).

The next theorem shows that the steplength α_k required to obtain the iterate ^(k+1) can be found
after a finite number of steps and that it is bounded away from zero.


Let δ∈ (0,1). There exist α> 0 independent of k and an integer j_k ≥ 0
such that for α_k = δ^  j_k

    F(^(k)_α) ≤  F(^(k)) + ηα_k ∇ F (^(k))^T (^(k) - ^(k)),
              
          α_k ≥α.




For F has Lipschitz continuous gradient, by applying Lemma <ref> we get

    F(^(k)_α)  ≤ F (^(k)) + α_k ∇ F (^(k))^T (^(k) - ^(k)) +
       α_k^2 L /2^(k) - ^(k)^2 ,

where L is the Lipschitz constant of ∇ F. Then, (<ref>) holds if we find
α_k such that

    ∇ F (^(k))^T (^(k) - ^(k)) + α_k  L /2^(k) - ^(k)^2
       ≤η∇ F (^(k))^T (^(k) - ^(k)),

or, equivalently,

    (1 - η) ∇ F (^(k))^T (^(k) - ^(k)) +
        α_k L /2^(k) - ^(k)^2 ≤ 0 .

From ∇ F_k(^(k)) = ∇ F(^(k)), the strong convexity of F_k and
step <ref> of Algorithm <ref>, it follows that

    ∇ F(^(k))^T (^(k) - ^(k))
                =      ∇ F_k(^(k))^T (^(k) - ^(k)) 
       ≤    F_k(^(k)) - F_k(^(k)) - γ/2^(k) - ^(k)^2
               
       ≤    - γ/2^(k) - ^(k)^2 ,

where γ is the strong convexity parameter of F_k. Thus, (<ref>)
holds for any α_k such that

    γ/2 (η - 1)  ^(k) - ^(k)^2 +
       α_k L /2^(k) - ^(k)^2 ≤ 0 .





By choosing the first nonnegative integer j_k such that

    δ^  j_k≤min{ 1, γ  (1 - η)/L}

and setting

    α= min{ 1, δ γ  (1 - η)/L}


we get the thesis.


Now we prove that the sequence generated by ACQUIRE  has a subsequence converging to
a solution of problem (<ref>). Because of the convexity of F, it is sufficient to prove that the
subsequence converges to a constrained stationary point of F.


Let {^(k)} be the sequence generated by Algorithm <ref>. Then there exists a subsequence
{^(k_j)} such that

    lim_k_j →∞^(k_j) = ,

where ∈ is such that ∇_ F() = 0.
Furthermore, any limit point  of {^(k)} is such that ∇_ F() = 0.


Let α_k = δ^  j^k, where j_k is given in Theorem <ref>.
By (<ref>) and (<ref>) we have

    F(^(k+1)) - F(^(k)) ≤ -  α_k η γ/2 ^(k) - ^(k)^2 ≤ 0 ;

then { F(^(k)) } is convergent, 
and the coercivity of F implies that {^(k)} is bounded. Since α_k ≥α> 0, we have that

    lim_k →∞^(k) - ^(k) = 0

and {^(k)} is bounded. This, together with ^(k) - ^(k)≤^(k)  - ^(k) + ^(k) - ^(k) and
the first inequality in (<ref>), implies that

    lim_k →∞^(k) - ^(k) = 0

and hence {^(k)} is bounded. Passing to subsequences, we have

    lim_k_j →∞^(k_j) = lim_k_j →∞^(k_j) = ∈.

Since the projection onto a nonempty closed convex set is nonexpansive, we get

    ∇_ F (^(k_j))  =
       ∇_ F (^(k_j)) - ∇_ F_k_j (^(k_j)) ≤∇ F (^(k_j)) - ∇ F_k_j (^(k_j))  ,

and,  by using (<ref>),

    lim_k_j →∞∇_ F (^(k_j))  =
        lim_k_j →∞∇ F (^(k_j)) - ∇ F_k_j (^(k_j))  = 0.

Then, for the lower semicontinuity of ∇_ F, we have

    ∇_ F() = 0̱.

If  is any limit point of {^(k)}, then ∈
and, by exploiting (<ref>) and passing to subsequences, we have

    lim_k_r →∞^(k_r) = lim_k_r →∞^(k_r) = ∈.

By reasoning as above we get

    ∇_ F() = 0̱,

which concludes the proof.


We note that ACQUIRE fits into the very general algorithmic framework presented
in <cit.> and hence Theorem <ref>
could be derived by specializing and adapting the convergence theory of that framework, taking into
account the specific properties of the functions D_KL() and TV_μ () and their quadratic approximations
D_KL^(k) () and TV_μ^(k) (), and the line search used. However, for the sake
of clarity and self-consistency, we decided to prove the convergence of Algorithm <ref> from scratch.

Now we show that if the objective function is strictly convex, the whole sequence {^(k)}
converges to the minimizer of problem (<ref>).


Assume that the function F is strictly convex. Then the sequence {^(k)} generated
by Algorithm <ref> converges to a point ∈ such that ∇_ F() = 0.


We follow the line of the proof of Lemma 2 in <cit.>.
By Theorem <ref> we know that there exists a limit point  of {^(k)}
such that  ∇_ F() = 0. Since F is strictly convex,  is the optimal solution
of problem (<ref>). We must prove that {^(k)} converges to .

From α_k ≤ 1 it follows that ^(k+1) - ^(k)≤^(k) - ^(k) and, by (<ref>),

    lim_k →∞^(k+1) - ^(k) = 0.

Since  is a strict minimizer, there exists δ > 0 such that F() < F() for all ∈
such that 0 <  - ≤δ. For all ε∈ (0, δ), it follows from Theorem <ref> that
the set B = {∈ : δ≤ - ≤ε} does not contain any limit point of {^(k)};

thus, there exists k_0 such that ^(k)∉B for all k > k_0 .

Let k_1 ≥ k_0 such that, for all k > k_1,

    ^(k+1) - ^(k)< δ - ε.

Let K be the set of indices defining a subsequence of {^(k)} converging to . There exists k ∈ K, k > k_1, such that

    ^(k) -  < ε ,

and hence

    ^(k+1) - ≤^(k+1) - ^(k) + ^(k) -  < δ - ε + ε = δ.

Since ^(k+1)∉B, we get

    ^(k+1) -< ε.

By the same argument we can prove that ^(k+j) -< ε implies ^(k+j+1) -< ε,
and hence, by induction, we have

    ^(k+j) -< ε   j.

Since ε is arbitrary, the thesis holds.


We conclude this section by showing that the stopping criterion (<ref>)
can be used to determine ^(k) at step <ref> of ACQUIRE.


Assume that (<ref>) holds for some θ∈ (0,1). Then, there exists {ε_k },
with ε_k>0 and lim_k →∞ε_k = 0, such that (<ref>) holds.



First we recall that - ∇ F_k() = ∇_ F_k() + P_N_ () (- ∇ F_k()) (see (<ref>)).
Since F_k is strongly convex with parameter γ and ^(k) is the
solution of problem (<ref>), we have

    γ/2^(k) -^(k)^2    ≤   
                      ( ∇ F_k(^(k)) - ∇ F_k(^(k)) )^T
                      ( ^(k) - ^(k) ) 
    
            =     ( ∇_ F_k(^(k)) )^T  ( ^(k) - ^(k) ) +
                  P_N_ (^(k)) ( - ∇ F_k(^(k)) )^T  ( ^(k) - ^(k) ) 
        +     P_N_ (^(k)) ( - ∇ F_k(^(k)) )^T  ( ^(k) - ^(k) ).

Since ^(k) - ^(k) belongs to the tangent cone at ^(k) and
^(k) - ^(k) belongs to the tangent cone at ^(k),

we get

    γ/2^(k) -^(k)^2 ≤
        ( ∇_ F_k(^(k)) )^T  ( ^(k) - ^(k) ) ≤∇_ F_k(^(k)) ^(k) - ^(k) .

It follows that

    ^(k) -^(k)≤2/γ∇_ F_k(^(k))   ;

thus, by requiring that

    ∇_ F_k(^(k))  ≤θ^k ∇_ F_k(^(0))

and setting ε_k = θ^k (2/γ) ∇_ F_k(^(0)), we get

    ^(k) -^(k)≤ε_k.












§ NUMERICAL EXPERIMENTS


ACQUIRE was implemented in MATLAB, using as inner solver the scaled gradient projection (SGP)
method proposed in <cit.>, widely applied in the solution of image restoration problems.
In particular, the implementation of SGP provided by the  MATLAB code,
available from <http://www.unife.it/prin/software>, was exploited.

The SGP iteration applied to problem (<ref>) reads:

    ^(j+1) = ^(j) + ρ_j ( P_,   C_j^-1( ^(j) - ν_j C_j ∇ F_k (^(j)) ) - ^(j)) ,

where ^(0) = ^(k), ρ_j is a line-search parameter ensuring that ^(j+1) satisfies a sufficient decrease
condition, ν_j is a suitably chosen steplength, C_j is a diagonal positive definite matrix with diagonal entries bounded
independently of j, and P_,   C_j^-1 is the projection operator onto  with respect to the norm induced
by the matrix C_j^-1 (the dependence on k has been neglected for simplicity).
Several efficient rules can be exploited to define the steplength ν_j for the quadratic problem (<ref>)  (see, e.g.,
<cit.> and the references therein).
In particular, SGP uses a modification of the ABB_ min adaptive Barzilai-Borwein steplength defined
in <cit.>, which takes into account the scaling matrix C_j (see <cit.>
for details); according to the analysis in <cit.>, this steplength appears very effective.
Since the steplength is computed by taking into account a certain number, say q, of suitable previous steplengths,
we modified  to avoid resetting the steplength each time the code was called, and to compute it by using q steplengths
from the previous call. The diagonal scaling matrix C_j was set as in <cit.> and
q was chosen equal to its defaul value in , i.e., q = 3.
The SGP iterations were stopped according to (<ref>). For all the tests considered here, we found experimentally
that θ = 0.1 worked well in the first iterations of ACQUIRE; on the other hand, criterion (<ref>) with this value
of θ soon becomes demanding, and fixing also a maximum number inner iterations was a natural choice. Setting this number
to 10 was effective in our experiments. Defaults were used for the remaining features of .

















The parameter γ in (<ref>) was set equal to 10^-5. The nonmonotone line search proposed
in <cit.> was implemented at line <ref>
of Algorithm <ref>, with memory length equal to 5, η = 10^-5, and δ = 0.5.
ACQUIRE was stopped using the following criterion

    ^(k+1) - ^(k)≤Tol ^(k) ,

i.e., when the relative change in the restored image went below a certain threshold.

ACQUIRE was compared with five state-of-the-art methods: PDAL, SGP,  SPIRAL-TV, SPLIT BREGMAN and VMILA.
By PDAL we denote our MATLAB implementation of the primal-dual algorithm proposed in <cit.>,
where we replaced the Chambolle-Pock algorithm <cit.> by the more efficient Primal Dual Algorithm with Linesearch introduced
in <cit.>. Concerning the parameters of PDAL, following <cit.> we set μ=0.7,
δ=0.99 and  β = 25. The initial steplength was chosen as τ = √(2/ω), where ω is an underestimate of M^T M

and M = [ A^T   D_1^T  …  D_n^T ]^T is the matrix linking the primal and dual variables.
SPIRAL-TV is the proximal-gradient method presented in <cit.>; a MATLAB
implementation of it is available from <http://drz.ac/code/>.
By SPLIT BREGMAN we denote a version of the method proposed in <cit.>,
which was specialized for problem (<ref>) <cit.>
and implemented in the MATLAB code  available from <http://dev.ipol.im/ getreuer/code/>.
Finally, VMILA is the variable-metric inexact line-search proximal-gradient method described in <cit.>,
whose MATLAB implementation can be found at
<http://www.oasis.unimore.it/site/home/software.html>. In all the methods, the stopping criterion (<ref>)
was applied. SGP was run with the same setting of parameters used to solve the subproblems in ACQUIRE.
For SPIRAL-TV, SPLIT BREGMAN and VMILA, the default values of the parameters were used.

PDAL, SPIRAL-TV, SPLIT BREGMAN and VMILA do not require any smooth approximation of TV and were
run directly on problem (<ref>). Therefore, our comparison also provides some insight into
the effects of using a smoothed version of TV. ACQUIRE was run with and without the flux constraint, i.e., using
both feasible sets _1 and _2 – see (<ref>) and (<ref>). However, since
the use of the flux constraint did not lead to any significant difference
in the restored images, and this constraint was not available in the implementations of  SPIRAL-TV, SPLIT BREGMAN
and VMILA, we report only the results for = _1.

As already noted, when the matrix A represents a convolution, the matrix-vector products involving the matrices A
and A^T can be performed by using fast algorithms. This is the case for all the experiments considered in this work.
Since periodic boundary conditions were considered for all the images used as test problems, the matrix-vector products were performed
by exploiting the MATLAB FFT functions  and .

In order to build the test problems used in the experiments, four reference images were chosen:
cameraman, micro, phantom and satellite,
shown in Figure <ref>. The cameraman image,
available in the MATLAB Image Processing Toolbox,
is widely used in the literature since it contains both sharp edges and flat regions
and presents a nice mixture of smooth and nonsmooth regions;

micro is the confocal microscopy phantom described in <cit.>;
phantom is the famous Shepp-Logan brain phantom described in <cit.>;
finally, the satellite image comes from the RestoreTools package <cit.>.
The size of cameraman, phantom and satellite is 256 × 256, while the size of
micro is 128 × 128.

A first set of test problems, T1, was obtained by convolving each reference image with a Gaussian PSF
and corrupting the resulting image with Poisson noise. A further set of test problems, T2, was built by
convolving some of the images with a motion blur PSF and an out-of-focus PSF, and then introducing Poisson noise.
Details about the PSFs and the Poisson noise are given in Subsections <ref> and <ref>,
where the results of numerical experiments performed by using the corrupted images are also reported.

All the experiments were carried out on a 2.5 GHz Intel Core i7 processors with 16 GB of RAM,
4 MB of L3 cache and the macOS 10.13.6 operating system, using MATLAB R2018b.




 §.§ Results on images with Gaussian blur


The Gaussian blur PSF for constructing the test set T1 was computed by using
the function  from <cit.>,
choosing the variance σ as specified in Table <ref>.
In order to take into account the existence of some background emission,
10^-10 was added to all the pixels of the blurred image; obviously, the vector  in D_KL()
was set as = 10^-10. The Poisson noise was introduced with the
function  from the MATLAB Image Processing Toolbox. Note that for this type of noise,
which affects the photon counting process, the Signal-to-Noise Ratio (SNR) is usually estimated by

    SNR = 10log_10(N_exact/√(N_exact + N_background)),

where N_exact and N_background are the total number of photons in the exact image to be recovered
and in the background term, respectively. Therefore, in order to obtain noisy and blurred images with SNR equal to 35 and 40,
the intensities of the reference images were suitably pre-scaled. The resulting images are shown
in Figures <ref>-<ref> (left columns).

The regularization parameter λ was set by trial and error, as described next.
The search for a suitable value of λ was carried out separately
for the minimization problem (<ref>), which uses the original TV, and the minimization
problem (<ref>), which uses the smoothed TV. In the latter case, ACQUIRE was run
several times on each test image, for 25 seconds and with Tol = 0 (see (<ref>)),
slowly varying the value of λ at each execution. The value of λ corresponding
to the smallest relative error at the last iteration, was chosen to perform the experiments with ACQUIRE
and SGP.  Note that, by running SGP for more than 25 seconds, we also verified that the selected value
of λ was suited to SGP too. The relative error was
computed as ^(k) - ^*  / ^*, where ^* denotes the original image.
The values of λ for problem (<ref>)
were set using the same strategy applied to (<ref>). In this case, instead of ACQUIRE,
for each image we used the method that appeared more efficient among PDAL, SPIRAL-TV, SPLIT-BREGMAN
and VMILA, on the basis of preliminary experiments. All the values of λ are reported in Table <ref>.
The same value of the regularization parameter was determined for both the original and the smoothed TV,
except for the satellite image; however, very close values of λ were obtained in this case.

The parameter μ in the smoothed version of TV was set as μ = 10^-2, which, by numerical experiments,
appeared to achieve a good tradeoff between approximation accuracy and computational effort, for all the test problems.
Indeed, as the value of μ decreases, TV_μ becomes a more accurate approximation of TV, but the
condition number of its Hessian increases. For both versions of TV,  each corrupted image was scaled by division
by its largest intensity value before applying any method; the scaled image was also used as starting guess, i.e., we set
^(0) =. We also performed experiments by setting ^(0) as the vector with entries equal to the
flux of the scaled image divided by the number of pixels of the image, but we could not see any significant difference in the results.

ACQUIRE was compared with all the other methods on the test problems previously described,
in terms of accuracy and execution time.
Six values of Tol were considered, Tol = 10^-2, 10^-3, …, 10^-7,
with the aim of assessing the behavior of ACQUIRE with different accuracy requirements
and getting useful information for the effective use of an automatic stopping rule.
A maximum execution time of 25 seconds was also set for all the methods.

Figures <ref> and <ref> show the relative errors and the execution times of each method,
in seconds, versus the stopping tolerances, for the problems with SNR = 35 and SNR = 40, respectively.
The images obtained with ACQUIRE and corresponding to the smallest
errors are shown in Figures <ref>-<ref> (right columns).
Further details concerning all the methods
are given in Tables <ref> and <ref>, where we report the
smallest errors, the iterations performed to achieve them, the corresponding execution times and tolerances,
and the values of MSSIM for the restored images corresponding to the smallest errors.  MSSIM is
a structural similarity measure index <cit.>
which is related to the perceived visual quality of the image; the higher its value, the better
the perceived similarity between the restored and original images.

We see that ACQUIRE generally does not need small tolerances to achieve small errors, because of its fast progress in the first
iterations, which produces large changes in the iterate. We note that in four test cases it reaches its minimum error
with Tol = 10^-3; this is consistent with the exploitation of second-order information to build the quadratic
model at each iteration. SGP generally achieves errors comparable with those of ACQUIRE, but its progress at each
iteration is slower, and hence it often requires smaller tolerances to avoid stopping prematurely. On the other hand,
a single iteration of ACQUIRE requires more time than an iteration of SGP, and the former method may be either
faster or slower than the latter in obtaining small errors. PDAL is able to achieve errors comparable with those of
ACQUIRE, but it generally requires smaller tolerances and larger times. VMILA is very efficient on both
instances of the cameraman problem and on the phantom problem with SNR =35, where it is faster than ACQUIRE
or comparable with it. However, there are some problems where VMILA makes very little progress in the first
iterations, leading to very premature stops, as shown by the almost constant execution times in the pictures.
The remaining methods are generally less efficient than the previous ones, because of their very slow progress in
reducing the error. We note that the errors of ACQUIRE and SGP show a light semiconvergence for some problems.
We were not able to completely remove this behavior by increasing the regularization parameter without significantly
deteriorating the visual quality of the image and decided to keep the value of λ determined by the
procedure previously described. Finally, we note that the values of MSSIM corresponding to ACQUIRE
confirm that in most cases this method is able to provide better or similar quality images in comparison with
the other methods.




 §.§ Results on images with moving and out-of-focus blurs


In order to understand if the previous behavior of ACQUIRE also holds for blurs different from the Gaussian
one, and to further compare ACQUIRE with the other methods, we built the test set T2. We introduced in two reference images,
cameraman and satellite, the motion blur and and the out-of-focus blur, which simulate the linear motion of a camera
and the out-of-focus effect, respectively. Both blurs were computed by using the Matlab function . 
Poisson noise with SNR equal to 35 and 40 was introduced in the blurred images, as in the case of Gaussian blur.
The length and the angle (in degrees) of the motion,  and φ, and the radius of the disk kernel
for the out-of-focus effect, , are specified in Table <ref>. The values of λ,
obtained with the procedure described in Subsection <ref>, are reported in the same table.
In this case, small differences can be observed between the values of the regularization parameter
corresponding to the original and the smoothed TV. The parameter μ in the smoothed version of TV
was set again as μ = 10^-2, on the basis of numerical experiments. Each noisy and blurred image
was scaled as in the case of Gaussian blur. The vector with entries equal to the flux of the scaled image
divided by the number of pixels of the image was used as starting guess, because, with the motion
and out-of-focus blurs, this choice generally appeared more effective than the choice of the scaled image.

The error and time histories shown in Figures <ref> and <ref>
confirm the behavior of ACQUIRE: it is able to strongly reduce the error in the first iterations and hence,
although its single iteration is usually more expensive than a single iterations of the other methods, it
is competitive with those methods. Furthermore, ACQUIRE allows an effective use
of an automatic stopping rule, avoiding premature stopping, which may happen with other methods.
This is confirmed by the data in Tables <ref> and <ref>,
which report, for all the methods, the smallest errors and the corresponding MSSIM values, number of iterations,
execution times and tolerances. The images restored by ACQUIRE and corresponding to the smallest errors are shown
in the right columns of Figures <ref>-<ref>). 





§ CONCLUSIONS


We proposed ACQUIRE, a method for TV-based restoration of images corrupted by Poisson noise, modeled
by (<ref>). ACQUIRE is a line-search method which considers a smoothed version of TV and computes
the search directions by minimizing quadratic models built by exploiting second-order information about the
objective function, which is usually not taken into account in methods for problem (<ref>).
We proved that the sequence generated by our method has a subsequence converging to a minimizer of
the smoothed problem (<ref>) and that any limit point is a minimizer; furthermore, if the
problem is strictly convex, the whole sequence is convergent.
We note that convergence holds without requiring the exact minimization of the quadratic models;
low accuracy in this minimization can be used in practice, as shown by the numerical results.

Computational experiments on reference test cases showed that the exploitation of second-order information
is beneficial, since it generally leads to a significant reduction of the reconstruction error in the first iterations,
Furthermore, the capability of achieving a tradeoff between accuracy and efficiency makes ACQUIRE
competitive with well-established methods for TV-based Poisson image restoration.

Finally, we observe that our approach can be extended to other regularization models,
such as high-order Total Variation <cit.> and Total Generalized Variation <cit.>, which have been proposed to reduce the staircase effect of TV and
retain the fine details of the image.




§ ACKNOWLEDGMENTS

We wish to thank the anonymous reviewers for their insightful comments and useful suggestions, which helped
us improve the quality of our work.


99
	
	
	BarnardBilheuxToops2018
	R.-C. Barnard, H. Bilheux, and T. e. a. Toops, Total variation-based
		neutron computed tomography, Review of Scientific Instruments, 89 (2018), 053704.
	
	BarzilaiBorwein1988
	J. Barzilai and J. M. Borwein, Two-point step size gradient
		methods, IMA Journal of Numerical Analysis, 8 (1988), pp. 141–148.
	
	BeckTeboulle2009
	A. Beck and M. Teboulle, Fast gradient-based algorithms for
		constrained total variation image denoising and deblurring problems, IEEE
	Transactions on Image Processing, 18 (2009), pp. 2419–2434.
	
	BerteroBoccacciDesideraVicidomini2009
	M. Bertero, P. Boccacci, G. Desiderà, and G. Vicidomini, Image
		deblurring with Poisson data: from cells to galaxies, Inverse Problems, 25
	(2009), 123006.
	
	BerteroBoccacciTalentiZanellaZanni2010
	M. Bertero, P. Boccacci, G. Talenti, R. Zanella, and L. Zanni, A
		discrepancy principle for poisson data, Inverse Problems, 26 (2010), 105004.
	
	BerteroLanteriZanni2008
	M. Bertero, H. Lantéri, and L. Zanni,
	Iterative image reconstruction: a point of view,
	in “Mathematical Methods in Biomedical Imaging and Intensity-Modulated Radiation Therapy (IMRT)”,
	pp. 37–63 (Y. Censor, M. Jiang and A. K. Louis eds.), CRM Series, Edizioni della Normale, 2008.
	
	Bertsekas1999
	D. P. Bertsekas, Nonlinear Programming, Athena Scientific, Belmont,
	MA, USA, 1999.
	
	BirginKrejicMartinez2003
	E. G. Birgin, N. Krejić, and J. M. Martínez, Globally convergent
		inexact quasi-Newton methods for solving nonlinear systems, Numerical
	Algorithms, 32 (2003), pp. 249–260.
	
	BonettiniLorisPortaPrato2016
	S. Bonettini, I. Loris, F. Porta, and M. Prato, Variable metric
		inexact line-search-based methods for nonsmooth optimization, SIAM Journal
	on Optimization, 26 (2016), pp. 891–921.
	
	BonettiniRuggiero2011
	S. Bonettini and V. Ruggiero, An alternating extragradient method
		for total variation-based image restoration from Poisson data, Inverse
	Problems, 27 (2011), 095001.
	
	BonettiniZanellaZanni2009
	S. Bonettini, R. Zanella, and L. Zanni, A scaled gradient projection
		method for constrained image deblurring, Inverse Problems, 25 (2009), 015002.
	
	BrediesHoller2014
	K. Bredies and M. Holler, Regularization of linear inverse problems
		with total generalized variation, Journal of Inverse and Ill-posed Problems,
	22 (2014), pp. 871–913.
	
	BrediesKunischPock2010
	K. Bredies, K. Kunisch, and T. Pock, Total generalized variation,
	SIAM Journal on Imaging Sciences, 3 (2010), pp. 492–526.
	
	CalamaiMore1987
	P. H. Calamai and J. J. Moré, Projected gradient methods for
		linearly constrained problems, Mathematical Programming, 39 (1987),
	pp. 93–116.
	
	CalamaiMore1987a
	P. H. Calamai and J. J. Moré, Quasi-Newton updates with
		bounds, SIAM Journal on Numerical Analysis, 24 (1987), pp. 1434–1441.
	
	Chambolle2004
	A. Chambolle, An algorithm for total variation minimization and
		applications, Journal of Mathematical Imaging and Vision, 20 (2004),
	pp. 89–97.
	
	ChambollePock2011
	A. Chambolle and T. Pock, A first-order primal-dual algorithm for
		convex problems with applications to imaging, Journal of Mathematical
	Imaging and Vision, 40 (2011), pp. 120–145.
	
	Condat2016
	L. Condat, Fast projection onto the simplex and the l_1 ball,
	Mathematical Programming, Series A, 158 (2016), pp. 575–585.
	
	DaiFletcher2006
	Y.-H. Dai and R. Fletcher, New algorithms for singly linearly
		constrained quadratic programs subject to lower and upper bounds,
	Mathematical Programming, Series A, 106 (2006), pp. 403–421.
	
	DeAsmundisdiSerafinoHagerToraldoZhang2014
	R. De Asmundis, D. di Serafino, W. W. Hager, G. Toraldo, and H. Zhang,
	An efficient gradient method using the Yuan steplength, Computational
	Optimization and Applications, 59 (2014), pp. 541–563.
	
	DeAsmundisdiSerafinoLandi2016
	R. De Asmundis, D. Di Serafino, and G. Landi, On the regularizing
		behavior of the SDA and SDC gradient methods in the solution of linear
		ill-posed problems, Journal of Computational and Applied Mathematics, 302
	(2016), pp. 81–93.
	
	DeAsmundisdiSerafinoRiccioToraldo2013
	R. De Asmundis, D. di Serafino, F. Riccio, and G. Toraldo, On
		spectral properties of steepest descent methods, IMA Journal of Numerical
	Analysis, 33 (2013), pp. 1416–1435.
	
	diSerafinoRuggieroToraldoZanni2018
	D. di Serafino, V. Ruggiero, G. Toraldo, and L. Zanni, On the
		steplength selection in gradient methods for unconstrained optimization,
	Applied Mathematics and Computation, 318 (2018), pp. 176–195.
	
	diSerafinoToraldoViolaBarlow2018
	D. di Serafino, G. Toraldo, M. Viola, and J. Barlow, A two-phase
		gradient method for quadratic programming problems with a single linear
		constraint and bounds on the variables, SIAM Journal on
	Optimization,  28 (2018), pp.  2809–2838.
	
	FacchineiLamparielloScutari2017
	F. Facchinei, L. Lampariello, and G. Scutari, Feasible methods for
		nonconvex nonsmooth problems with applications in green communications,
	Mathematical Programming, 164 (2017), pp. 55–90.
	
	FigueiredoBioucasDias2010
	M. A. T. Figueiredo and J. M. Bioucas-Dias, Restoration of
		Poissonian images using alternating direction optimization, IEEE
	Transactions on Image Processing, 19 (2010), pp. 3133–3145.
	
	Fletcher2012
	R. Fletcher, A limited memory steepest descent method, Mathematical
	Programming, Series A, 135 (2012), pp. 413–436.
	
	FrassoldatiZanniZanghirati2008
	G. Frassoldati, L. Zanni, and G. Zanghirati, New adaptive stepsize
		selections in gradient methods, Journal of Industrial and Management
	Optimization, 4 (2008), pp. 299–312.
	
	Getreuer2012
	P. Getreuer, Total variation deconvolution using split Bregman,
	Image Processing On Line, 2 (2012), pp. 158–174.
	
	GoldsteinOsher2009
	T. Goldstein and S. Osher, The split Bregman method for
		L1-regularized problems, SIAM Journal on Imaging Sciences, 2 (2009),
	pp. 323–343.
	
	GrippoLamparielloLucidi1986
	L. Grippo, F. Lampariello, and S. Lucidi, A nonmonotone line search
		technique for Newton's method, SIAM Journal on Numerical Analysis, 23
	(1986), pp. 707–716.
	
	HarmanyMarciaWillet2012
	Z. T. Harmany, R. F. Marcia, and R. M. Willett, This is
		SPIRAL–TAP: Sparse Poisson Intensity Reconstruction Algorithms
		– Theory and Practice, IEEE Transactions on Image Processing, 21
	(2012), pp. 1084–1096.
	
	Herman2009
	G. Herman, Fundamentals of Computerized Tomography, Springer,
	Berlin, 2009.
	
	JohnssonHuangChan1998
	E. Jonsson, S.-C. Huang, and T. Chan, Total Variation regularization
		in positron emission tomography, Tech. Report 98–48, UCLA, CAM-Report, 1998.
	
	LiMalgouyresZeng2017
	Z. Li, F. Malgouyres, and T. Zeng,
	Regularized non-local Total Variation and application in image restoration
	Journal of Mathematical Imaging and Vision, 59 (2017), pp. 296–317.
	
	LiuHuangLvWang2017
	J. Liu, T.-Z. Huang, X.-G. Lv, and S. Wang, High-order total
		variation-based Poissonian image deconvolution with spatially adapted
		regularization parameter, Applied Mathematical Modelling, 45 (2017),
	pp. 516–529.
	
	PiccolominiColiMorottiZanni2017
	E. Loli Piccolomini, V. L. Coli, E. Morotti, and L. Zanni, Reconstruction of 3D X-ray CT images from reduced sampling by a scaled
		gradient projection algorithm, Computational Optimization and Applications, 71 (2018),
		pp. 171–191. 
	
	MalitskyPock2018
	Y. Malitsky and T. Pock, A first-order primal-dual algorithm with
		linesearch, SIAM Journal on Optimization, 28 (2018), pp. 411–432.
	
	MoreToraldo1991
	J. J. Moré and G. Toraldo, On the solution of large quadratic
		programming problems with bound constraints, SIAM Journal on Optimization, 1
	(1991), pp. 93–113.
	
	MotaMatelaOliveraAlmeida2016
	A.-M. Mota, N. Matela, N. Olivera, and P. Almeida, 3D total
		variation minimization filter for breast tomosynthesis imaging, Proc. of
	13th International Workshop IWDM,  (2016), pp. 484–492.
	
	NagyPalmerPerrone2004
	J. G. Nagy, K. Palmer, and L. Perrone, Iterative methods for image
		deblurring: a Matlab object-oriented approach, Numerical Algorithms, 36
	(2004), pp. 73–93.
	
	PaninZengGullberg1999
	V. Y. Panin, G. L. Zeng, and G. T. Gullberg, Total variation
		regulated EM algorithm, in 1998 IEEE Nuclear Science Symposium Conference
	Record, 1998, p. 6359010.
	
	PapafitsorosSchonlieb2014
	K. Papafitsoros and C.-B. Schönlieb, A combined first and second
		order variational approach for image reconstruction, Journal of Mathematical
	Imaging and Vision, 48 (2014), pp. 308–338.
	
	Pawley1996
	J. B. Pawley, Handbook of Biological Confocal Microscopy, Plenum
	Press, New York, 2nd ed., 1996.
	
	RodriguezWohlberg2009
	P. Rodríguez and B. Wohlberg, Efficient minimization method for
		a generalized total variation functional, IEEE Transactions on Image
	Processing, 18 (2009), pp. 322–332.
	
	RudinOsherFatemi1992
	L. I. Rudin, S. Osher, and E. Fatemi, Nonlinear total variation
		based noise removal algorithms, Physica D, 60 (1992), pp. 259–268.
	
	SarderNehorai2006
	P. Sarder and A. Nehorai, Deconvolution method for 3-D
		fluorescence microscopy images, IEEE Signal Processing Letters, 23 (2006),
	pp. 32–45.
	
	SawatzkyBruneWuebbelingKoestersSchaefersBurger2008
	A. Sawatzky, C. Brune, F. Wübbeling, T. Kösters, K. Schäfers, and
		M. Burger, Accurate EM-TV algorithm in PET with low SNR, in 2008
	IEEE Nuclear Science Symposium Conference Record, 2008,
	<https://doi.org/10.1109/NSSMIC.2008.4774392>.
	
	SetzerSteidlTeuber2010
	S. Setzer, G. Steidl, and T. Teuber, Deblurring Poissonian images
		by split Bregman techniques, Journal of Visual Communication and Image
	Representation, 21 (2010), pp. 193–199.
	
	SheppLogan1974
	L. A. Shepp and B. Logan, The Fourier reconstruction of a head
		section, IEEE Transactions on Nuclear Science, 21 (1974), pp. 21–43.
	
	SheppVardi1982
	L. A. Shepp and Y. Vardi, Maximum likelihood reconstruction for
		emission tomography, IEEE Transactions on Medical Imaging, 1 (1982),
	pp. 113–122.
	
	StarckMurtagh2006
	J.-L. Starck and F. Murtagh, Astronomical Image and Data Analysis,
	Springer, 2006.
	
	VogelOman1996
	C. R. Vogel and M. E. Oman, Iterative methods for total variation
		denoising, SIAM Journal on Scientific Computing, 17 (1996), pp. 227–238.
	
	WangBovikSheikhSimoncelli2004
	Z. Wang, A. C. Bovik, H. R. Sheikh, and E. P. Simoncelli, Image
		quality assessment: from error visibility to structural similarity, IEEE
	Transactions on Image Processing, 13 (2004), pp. 600–612.
	
	WeissBlancFeraudAubert2009
	P. Weiss, L. Blanc-Féraud, and G. Aubert, Efficient schemes for
		total variation minimization under constraints in image processing, SIAM
	Journal on Scientific Computing, 31 (2009), pp. 2047–2080.
	
	WenChanZeng2016
	Y.-W. Wen, R. H. Chan, and T.-Y. Zeng, Primal-dual algorithms for
		total variation based image restoration under poisson noise, Science China
	Mathematics, 59 (2016), pp. 141–160.
	
	WillettNowak2003
	R. M. Willett and R. D. Nowak, Platelets: a multiscale approach for
		recovering edges and surfaces in photon limited medical imaging, IEEE
	Transactions on Medical Imaging, 22 (2003), pp. 332–350.
	
	ZanellaBoccacciZanniBertero2009
	R. Zanella, P. Boccacci, L. Zanni, and M. Bertero, Efficient
		gradient projection methods for edge-preserving removal of Poisson noise,
	Inverse Problems, 25 (2009), 045010.
	
	Zarantonello1971
	E. H. Zarantonello, Projections on convex sets in Hilbert space
		and spectral theory, in Contributions to Nonlinear Fanctional Analysis,
	E. H. Zarantonello, ed., Academic Press, New York, NY, USA, 1971.
	
	ZhangHuNagy2018
	J. Zhang, Y. Hu, and J. Nagy, A scaled gradient method for digital
		tomographic image reconstruction, Inverse Problems and Imaging, 12 (2018),
	pp. 239–259.

	









































































