Chao Wang
email: chaowang.hk@gmail.com

Raymond Chan
email: rchan@math.cuhk.edu.hk

Mila Nikolova
email: nikolova@cmla.ens-cachan.fr

Robert Plemmons
email: plemmons@wfu.edu

Sudhakar Prasad
email: sprasad@unm.edu

Nonconvex optimization for 3D point source localization using a rotating point spread function *

Keywords: nonconvex optimization algorithms, 3D localization, space debris, point spread function, image rotation, image processing AMS subject classifications. 65K10, 65F22, 90C26

We consider the high-resolution imaging problem of 3D point source image recovery from 2D data using a method based on point spread function (PSF) engineering. The method involves a new technique, recently proposed by S. Prasad, based on the use of a rotating PSF with a single lobe to obtain depth from defocus. The amount of rotation of the PSF encodes the depth position of the point source. Applications include high-resolution single molecule localization microscopy as well as the problem addressed in this paper on localization of space debris using a space-based telescope. The localization problem is discretized on a cubical lattice where the coordinates of nonzero entries represent the 3D locations and the values of these entries the fluxes of the point sources. Finding the locations and fluxes of the point sources is a large-scale sparse 3D inverse problem. A new nonconvex regularization method with a data-fitting term based on Kullback-Leibler (KL) divergence is proposed for 3D localization for the Poisson noise model. In addition, we propose a new scheme of estimation of the source fluxes from the KL data-fitting term. Numerical experiments illustrate the efficiency and stability of the algorithms that are trained on a random subset of image data before being applied to other images. Our 3D localization algorithms can be readily applied to other kinds of depth-encoding PSFs as well.

1. Introduction. Imaging and localizing point sources with high accuracy in a 3D volume is an important and challenging task. An example is super-resolution 3D single-molecule localization [START_REF] Birk | Super-Resolution Microscopy: A Practical Guide[END_REF][START_REF] Von Diezmann | Three-dimensional localization of single molecules for super-resolution imaging and single-particle tracking[END_REF], which is an area of intense interest in biology (cell imaging, folding, membrane behavior, etc.), in chemistry (spectral diffusion, molecular distortions, etc.), and in physics (structures of materials, quantum optics, etc.). A different problem of interest in the space-surveillance community is that of localization of space debris using a space-based telescope. Since the optical wavelength is much shorter than the radio wavelength, optical detection and localization is expected to attain far greater precision than the more commonly employed radar systems. However, the shorter field depth of optical imaging systems may limit their performance to a shorter range of distances. An integrated system consisting of a radar system for performing radio detection, localization, and ranging of space debris at larger distances, which cues in an optical system when debris reach shorter distances, may ultimately provide optimal performance for detecting and tracking debris at distances ranging from tens of kilometers down to hundreds of meters.

A stand-alone optical system based on the use of a light-sheet illumination and scattering concept [START_REF] Englert | Optical orbital debris spotter[END_REF] for spotting debris within meters of a spacecraft has been proposed. A second system can localize all three coordinates of an unresolved, scattering debris [START_REF] Hampf | Optical technologies for the observation of low earth orbit objects[END_REF][START_REF] Wagner | Detection and laser ranging of orbital objects using optical methods[END_REF] by utilizing either parallex between two observatories or a pulsed laser ranging system or a hybrid system. For parallex, two observatories receive debris scattered optical signal simultaneously. For the pulsed laser, the ranging system is coupled to a single imaging observatory. The hybrid system utilizes both approaches in which the laser pulse transmitted from one of the two observatories is received at time-gated single-photon detectors with good parallax information at both the observatories. However, to the best of our knowledge there is no other proposal for a full 3D debris localization and tracking optical or optical-radar system working in the range of tens to hundreds of meters. Prasad [START_REF] Prasad | Innovations in space-object shape recovery and 3D space debris localization[END_REF] has proposed the use of an optical imager that exploits off-center image rotation to encode in a single image snapshot both the range z and transverse (x, y) coordinates of a swarm of unresolved sources such as small, sub-centimeter class space debris, which when actively illuminated can scatter a fraction of laser irradiance back into the imaging sensor.

Here we develop a promising non-convex optimization algorithm and compare its performance with other recent algorithms that can reconstruct the 3D positions and fluxes of a random collection of many point sources within the focal depth of a rotating-PSF-based imaging system from a single image dataset. We shall assume the data to be corrupted only by signal-dependent shot noise that is well described by a Poisson statistical model, as would be roughly characteristic of signal acquisition by an EM-CCD sensor in the photon counting (PC) mode [START_REF] Daigle | Photon counting with an EMCCD[END_REF].

1.1. Previous related work. The area of 3D object localization and imaging is getting increased attention in recent years. One method is to scan the 2D slices of information at different depths and then reconstruct the 3D image [START_REF] Prabhat | Simultaneous imaging of different focal planes in fluorescence microscopy for the study of cellular dynamics in three dimensions[END_REF]. Owing to the inefficiency and other limitations of standard PSFs, several optical modifications [START_REF] Birk | Super-Resolution Microscopy: A Practical Guide[END_REF][START_REF] Von Diezmann | Three-dimensional localization of single molecules for super-resolution imaging and single-particle tracking[END_REF] have been developed for these problems. There are three typical classes of methods: multifocus method, interferometric detection and point spread function (PSF) engineering. The multifocus method [START_REF] Juette | Three-dimensional sub-100 nm resolution fluorescence microscopy of thick samples[END_REF][START_REF] Prabhat | Simultaneous imaging of different focal planes in fluorescence microscopy for the study of cellular dynamics in three dimensions[END_REF] uses more than one focal plane simultaneously. The interferometric method [START_REF] Shtengel | Interferometric fluorescent super-resolution microscopy resolves 3D cellular ultrastructure[END_REF] extracts the locations of point sources from the interference fringes of a source signal coherently propagating in two opposed paths in the imaging instrument. PSF engineering [START_REF] Pavani | High-efficiency rotating point spread functions[END_REF][START_REF] Pavani | Three-dimensional, single-molecule fluorescence imaging beyond the diffraction limit by using a double-helix point spread function[END_REF][START_REF] Lew | Corkscrew point spread function for far-field three-dimensional nanoscale localization of pointlike objects[END_REF][START_REF] Prasad | Rotating point spread function via pupil-phase engineering[END_REF][START_REF] Shechtman | Optimal point spread function design for 3D imaging[END_REF][START_REF] Shechtman | Precise three-dimensional scan-free multiple-particle tracking over large axial ranges with tetrapod point spread functions[END_REF], which encodes the depth coordinate of point sources in a 3D scene into a single 2D snapshot, is based on choosing a phase pattern that makes the defocused image of a point source depthdependent without blurring it excessively. For example, this can be achieved by inserting a cylindrical lens into an optical system to get approximately elliptical Gaussian, astigmatic images whose axial orientation and size contain unique information about the source locations relative to the plane of best focus [START_REF] Huang | Three-dimensional super-resolution imaging by stochastic optical reconstruction microscopy[END_REF]. A more general way to create the requisite phase aberrations is to use phase masks. There are many different kinds of phase masks, which give rise to the double-helix PSF [START_REF] Pavani | Three-dimensional, single-molecule fluorescence imaging beyond the diffraction limit by using a double-helix point spread function[END_REF][START_REF] Pavani | High-efficiency rotating point spread functions[END_REF], the corkscrew PSF [START_REF] Lew | Corkscrew point spread function for far-field three-dimensional nanoscale localization of pointlike objects[END_REF], the Tetrapod PSFs [START_REF] Shechtman | Optimal point spread function design for 3D imaging[END_REF][START_REF] Shechtman | Precise three-dimensional scan-free multiple-particle tracking over large axial ranges with tetrapod point spread functions[END_REF] and the single-lobe PSF [START_REF] Prasad | Rotating point spread function via pupil-phase engineering[END_REF][START_REF] Kumar | PSF rotation with changing defocus and applications to 3D imaging for space situational awareness[END_REF]. Since the amount of rotation of the PSF encodes the depth information for the double-helix PSF, the corkscrew and the single-lobe PSF, we call them rotating PSFs. Tetrapod PSFs [START_REF] Shechtman | Optimal point spread function design for 3D imaging[END_REF][START_REF] Shechtman | Precise three-dimensional scan-free multiple-particle tracking over large axial ranges with tetrapod point spread functions[END_REF] represent another approach that is rather analogous to astigmatic imaging in that the source depth is encoded in the shape of the PSF. In this paper, we will only focus on Prasad's single-lobe rotating PSF [START_REF] Prasad | Rotating point spread function via pupil-phase engineering[END_REF]. Let us explain the movitivation for this choice.

Several algorithms have been developed for point source localization, for example in the area of single molecule localization in microscopy 3D DAOSTORM [START_REF] Babcock | A high-density 3D localization algorithm for stochastic optical reconstruction microscopy[END_REF] and 3D FALCON [START_REF] Min | 3D high-density localization microscopy using hybrid astigmatic/biplane imaging and sparse image reconstruction[END_REF]. However these two methods cannot be applied for localization using rotating PSFs, since they are based on fitting Gaussian functions, which works well only for astigmatism-based PSFs. Easy-DHPSF [START_REF] Lew | Easy-DHPSF open-source software for threedimensional localization of single molecules with precision beyond the optical diffraction limit[END_REF] is a popular approach for solving the double-helix PSF problem but it requires no overlap for PSF images which is often impossible to satisfy in reality. In [START_REF] Shuang | Generalized recovery algorithm for 3D super-resolution microscopy using rotating point spread functions[END_REF], the authors proposed a generalized method for various rotating PSFs by using 3D deconvolution with a regularization method. The optimization model is for an additive Poisson noise model that is only approximately correct even under highly limited brightness conditions. Owing to sparsity, a hybrid algorithm with matching pursuit and convex optimization was described in [START_REF] Barsic | Three-dimensional super-resolution and localization of dense clusters of single molecules[END_REF] without any theoretical outline. We therefore believe that it is important to propose an algorithm for the high-density case without the assumption of highly limited brightness condition. With this aim, we propose an optimization model for 3D localization and consider the single-lobe rotating PSF [START_REF] Prasad | Rotating point spread function via pupil-phase engineering[END_REF]. We emphasize, however, that our model is broadly applicable to a variety of PSFs.

Problem development.

In image formation modeling, the point spread function (PSF) is the imaging system's response to light from a point source. Our purpose is to extend the work of Prasad on PSF engineering [START_REF] Prasad | Rotating point spread function via pupil-phase engineering[END_REF], which proposed the use of depth-dependent image rotation, by developing a nonconvex optimization algorithm for 3D point source localization using such an imaging system. See Section 2 for details.

By imposing spiral phase retardation with a phase winding number that changes in regular integer steps from one annular zone to the next of an aperture-based phase mask, one can create an image of a point source that has an approximate rotational shape invariance, provided the zone radii have a square root dependence on their indices. Specifically, when the distance of the source from the aperture of such an imaging system changes, the off-center, shape-preserving PSF merely rotates by an amount roughly proportional to the source misfocus from the plane of best focus. The following general model based on the rotating PSF image describes the spatial distribution of image brightness for P point sources. The observed 2D image without noise is

(1.1) I 0 (x, y) = P i=1 H i (x -x i , y -y i )f i + b,
where b is the uniform background and H i (x -x i , y -y i ) is the rotating PSF for the i-th point source of flux f i and 3D position coordinates (x i , y i , z i ) with the depth information coded in each H i , and (x, y) is the position in image plane.

Here, we build a forward model for the problem based in part on the approach developed in [START_REF] Shuang | Generalized recovery algorithm for 3D super-resolution microscopy using rotating point spread functions[END_REF]. In order to estimate the 3D locations of the point sources, we assume their distribution is approximated by a discrete lattice X ∈ R m×n×d . The indices of the nonzero entries of X are the 3-dimensional locations of the point sources and the values at these entries correspond to the fluxes, i.e., the energy emitted by the illuminated point source. The 2D observed image G ∈ R m×n can be represented as

G = P (T (A * X ) + b) ,
where b is background signal and A * X is the convolution of X with the 3D PSF A. This 3D PSF is a cube which is constructed by a sequence of images with respect to different depths of the points. Each slice is the image corresponding to a point source at the origin in the (x, y) plane and at depth z. Here T is an operator for extracting the last slice of the cube A * X since the observed information is a snapshot, and P is the Poisson noise operator. Figure 1 shows the forward model for a specific PSF. In order to recover X , we need to solve a large-scale sparse 3D inverse problem given as follows:

min X X 0 s.t. D (X ) < (1.2)
where X 0 is the counting function which gives the number of nonzero entries in X and is also called the 0 norm. Here provides the noise level and D is a certain data-fitting term based on the noise model. This problem is NP-hard [START_REF] Natarajan | Sparse approximate solutions to linear systems[END_REF] and there are a number of algorithms and methods [START_REF] Foucart | A mathematical introduction to compressive sensing[END_REF] being developed to find sparse representation or a good approximate solution. Examples of such algorithms include greedy algorithms like orthogonal matching pursuit (OMP) [START_REF] Pati | Orthogonal matching pursuit: Recursive function approximation with applications to wavelet decomposition[END_REF], iterative hard thresholding algorithm (IHT) [START_REF] Blumensath | Iterative thresholding for sparse approximations[END_REF], which requires the knowledge of the exact number of sparse entries in the ground truth, and convex/nonconvex relaxation methods [START_REF] Nikolova | Efficient reconstruction of piecewise constant images using nonsmooth nonconvex minimization[END_REF][START_REF] Nikolova | Fast nonconvex nonsmooth minimization methods for image restoration and reconstruction[END_REF][START_REF] Nikolova | On 1 data fitting and concave regularization for image recovery[END_REF][START_REF] Xiao | On the convergence of nonconvex minimization methods for image recovery[END_REF][START_REF] Chen | Atomic decomposition by basis pursuit[END_REF][START_REF] Soubies | A continuous exact 0 penalty (CEL0) for least squares regularized problem[END_REF] where the 0 term is approximated and (1.2) is replaced by a regularization model. Here, we develop a convex/nonconvex relaxation approach to solving (1.2) for the Poisson noise model. In comparison to state-of-the-art methods, our approach is more efficient even for high density cases. We do not need to know the exact number of the ground truth point sources and, furthermore, we have observed less overfitting than other methods. That is to say that our estimated result recovers fewer false point sources.

1.3. Outline of the paper. The rest of the paper is organized as follows. In Section 2, we describe the physics model for Prasad's single lobe rotating PSF. In Section 3, we propose a nonconvex optimization model to solve the point source localization problem. In Section 4, our nonconvex optimization algorithm is developed, with a post-processing step for eliminating clustered false positive point sources. In addition, a new iterative scheme for estimating the flux values is also proposed in this section. Numerical experiments, including comparisons with other optimization models, are discussed in Section 5. Some concluding remarks are made in Section 6. 2. Physics model for single lobe rotating PSFs . Here, we consider a new technique, recently patented by Prasad in [START_REF] Kumar | PSF rotation with changing defocus and applications to 3D imaging for space situational awareness[END_REF][START_REF] Prasad | Rotating point spread function via pupil-phase engineering[END_REF], for applying a rotating PSF with a single lobe to obtain depth from defocus. Comparing to the traditional double-helix rotating PSF [START_REF] Pavani | Three-dimensional, single-molecule fluorescence imaging beyond the diffraction limit by using a double-helix point spread function[END_REF], this single lobe PSF uses photons from targets more efficiently and yields images that are less cluttered and confounding in high-density swarm scenarios than other techniques.

The amount of rotation of Prasad's PSF encodes the depth position of the point source. Denote A ζ as the PSF matrix for a point source with flux f = 1, transverse location r 0 = (x 0 , y 0 ) = 0, and defocus parameter ζ. Note that ζ is proportional to the object-distance δz from the in-focus object plane, according to the relation

(2.1) ζ = - πδzR 2 λl 0 (l 0 + δz) ,
where l 0 is the distance between the lens and the best focus point and R is the radius of the pupil. Here the imaging wavelength is denoted by λ. The two-dimensional image data matrix G ∈ R m×n with Poisson noise is generated as:

(2.2)

G = P P i=1 f i (A ζ i * δ(x i , y i )) + b ,
where b is the spatially uniform background, * is the convolution operation, and δ(x i , y i ) is the location matrix for the ith source, with elements,

(δ(x i , y i )) uv = 1, (u, v) = (x i , y i ), 0, otherwise.
The center of rotation of the image in the transverse plane, denoted by r I = (x, y), is related to the source transverse location, r 0 , in the object plane as

r I = - z I r 0 l 0 + δz ,
where z I is the distance between the image plane and the lens.

According to the Fourier optics model, the incoherent PSF for a clear aperture containing a phase mask with optical phase retardation, ψ(s), is given by

A ζ (s) = F -1 P (s)exp ι(ζs 2 -ψ(s)) 2 ,
where ι = √ -1, F -1 denotes the inverse Fourier transform, P (s) denotes the indicator function for the pupil of radius R, and s = (s, φ s ) is a scaled version of the image-plane position vector, r, namely s = r λz I /R . Here r is measured from the center of the geometric image point located at r I . For the single-lobe rotating PSF, ψ(s) is chosen to be the spiral phase profile defined as

ψ(s) = lφ s , for l -1 L ≤ s ≤ l L , l = 1, • • •, L,
in which L is the number of concentric annular zones in the phase mask.

Remark: For the discretized forward model (see Figure 1), we get the dictionary A by sampling depths at regular intervals in the range, ζ i ∈ [-πL, πL], over which the PSF performs one complete rotation about the geometric image center before it begins to break apart. The i-th slice of dictionary is denoted as

A ζ i .
3. Nonconvex optimization model . In general, for inverse problems the objective function includes a data-fitting term and a regularization term. The regularization term, representing prior information, is based on the desired properties of the solution, while the data-fitting term measures the discrepancy between the estimated and the observed images according to the noise model. Here we consider the Poisson noise model for which the data fitting term is the I-divergence. The model with a regularization term, µR(X ) is min

X ≥0 D KL (T (A * X ) + b 1, G) + µR(X ),
where the uniform background is denoted as b. Here we also use it to denote the matrix b 1, 1 ∈ R m×n , i.e., b 1 is the m × n matrix with all entries b. The data-fidelity term is the negative log-likelihood of the data. For Poisson noise, this term is the I-divergence which is also known as Kullback-Leibler (KL) divergence (see [START_REF] Le | A variational approach to reconstructing images corrupted by poisson noise[END_REF])

D KL (z, g) = g, ln g z + 1, z -g .
This is a large-scale sparse 3D inverse problem. For the regularization term, we consider a nonconvex nondifferentiable function (see [START_REF] Nikolova | Efficient reconstruction of piecewise constant images using nonsmooth nonconvex minimization[END_REF][START_REF] Nikolova | Fast nonconvex nonsmooth minimization methods for image restoration and reconstruction[END_REF][START_REF] Nikolova | On 1 data fitting and concave regularization for image recovery[END_REF][START_REF] Xiao | On the convergence of nonconvex minimization methods for image recovery[END_REF]), using specifically

R(X ) := m,n,d i,j,k=1 θ(X ijk ) = m,n,d i,j,k=1 |X ijk | a + |X ijk | ,
where a is fixed and determines the degree of nonconvexity, see Figure 2(a). Therefore, the minimization problem amounts to

(3.1) min X ≥0    1, T (A * X ) -G ln(T (A * X ) + b 1) + µ m,n,d i,j,k=1 |X ijk | a + |X ijk |    .
Here, θ(t) = lim →1 θ (t), where θ (t) = |t| a+ |t| . Since θ (t) represents the 1 norm when = 0, we see the process of increasing nonconvexity as increases from 0 to 1; see Figure 2(b). Moreover, we easily get the derivative of θ(t):

θ (t) = a sign(t) (a + |t|) 2
, when t = 0, and θ (0

+ ) = 1 a ,
where θ (0 + ) is the right derivative at zero. 4. Algorithm development. The minimization problem (3.1) involves nonconvex nondifferentiable optimization. Here we use an iterative reweighted 1 algorithm (IRL1) [START_REF] Ochs | On iteratively reweighted algorithms for nonsmooth nonconvex optimization in computer vision[END_REF] to solve this optimization problem. It is a majorization-minimization method which solves a series of convex optimization problem with a weighted-1 regularization term. It considers the problem (see Algorithm 3, in [START_REF] Ochs | On iteratively reweighted algorithms for nonsmooth nonconvex optimization in computer vision[END_REF])

min x∈X F (x) := F 1 (x) + F 2 (G(x)),
where X is the constraint set. F is a lower semicontinuous (lsc) function, extended, realvalued, proper, while F 1 is proper, lower-semicontinous, and convex and F 2 is coordinatewise nondecreasing, i.e. F 2 (x) ≤ F 2 (x + te i ) with x, x + te i ∈ G(X) and t > 0, where e i is the i-th canonical basis unit vector. The function F 2 is concave on G(X). The IRL1 iterative scheme

[25, Algorithm 3] is    w l = ∂F 2 (y), y = G(x l ), x l+1 = argmin x∈X F 1 (x) + w l , G(x) ,
where ∂ stands for subdifferential. For our problem (3.1), we can choose:

F 1 (X ) = 1, T (A * X ) -G log(T (A * X ) + b 1) ; F 2 (X ) = µ m,n,d i,j,k=1 X ijk a + X ijk ; G(X ) = |X |; X = {X | X ijk ≥ 0 for all i, j, k}.
Therefore, we compute the partial derivative of w l and get w l ijk = aµ (a+ X l ijk )

2 , ∀i, j, k. Here w l = 0 is finite, since a, µ = 0, and all X ijk ≥ 0 owning to the constraint X. According to [START_REF] Nikolova | Efficient reconstruction of piecewise constant images using nonsmooth nonconvex minimization[END_REF][START_REF] Harmany | This is spiral-tap: Sparse poisson intensity reconstruction algorithmstheory and practice[END_REF], these terms satisfy the requirements of the algorithm. Algorithm 4.1 gives the IRL1 for solving (3.1).

Algorithm 4.1 Iterative reweighted 1 algorithm (IRL1) for the rotating PSF problem

Input: X 0 ∈ R m×n×d and G ∈ R m×n . Set a and µ. Output: The solution X * which is the minimizer in the last outer iteration.

1: repeat

2:

Compute w l ijk = aµ (a+ X l ijk )

2 , ∀i, j, k;

3:
Given G, w l ijk , obtain X l by solving (4.1) X l = argmin

X ≥0    1, T (A * X ) -G log(T (A * X ) + b 1) + m,n,d i,j,k=1 w l ijk |X ijk |    . 4: until convergence 4.
1. Subproblem of IRL1. In Algorithm 4.1, we need to solve the subproblem (4.1), which can be solved by the alternating direction method of multipliers (ADMM) [START_REF] Boyd | Distributed optimization and statistical learning via the alternating direction method of multipliers[END_REF]. For this, we introduce two auxiliary variables, namely U 0 and U 1 , and transform (4.1) into min

U 1 ≥0,U 0 ,X    1, T U 0 -G log(T U 0 + b 1) + m,n,d i,j,k=1 w l ijk |(U 1 ) ijk | : U 0 = A * X , U 1 = X    .
The nonnegative constraint can be rewritten as an indicator function such that the objective function in this subproblem becomes

(4.2) 1, T U 0 -G log(T U 0 + b 1) + m,n,d i,j,k=1 w l ijk |(U 1 ) ijk | + 1 + (U 1 ),
where 1 + (X ) is the indicator function of the nonnegative constraint

1 + (X ) = 0, X ≥ 0, ∞, otherwise.
Therefore, the augmented Lagrangrian function L(U 0 , U 1 , X , η 0 , η 1 ) [START_REF] Nocedal | Numerical optimization[END_REF] for subproblem (4.1) is

L(U 0 , U 1 , X , η 0 , η 1 ) := 1, T U 0 -G log(T U 0 + b 1) + m,n,d i,j,k=1 w l ijk |(U 1 ) ijk | + β 0 2 U 0 -A * X -η 0 2 F + β 1 2 U 1 -X -η 1 2 F + 1 + (U 1 ),
where η 0 , η 1 ∈ R m×n×d are the Lagrange multipliers and β 0 , β 1 > 0. Here X F is the Frobenius norm of X , which is equal to the 2 norm of the vectorized X . Starting at X = X t , U 1 = U t 1 and U 2 = U t 2 , applying ADMM in [START_REF] Boyd | Distributed optimization and statistical learning via the alternating direction method of multipliers[END_REF] yields the iterative scheme

U t+1 0 = argmin U 0 1, T U 0 -G log(T U 0 + b 1) + β 0 2 U 0 -A * X t -η t 0 2 F (4.3a) U t+1 1 = argmin U 1 ≥0    m,n,d i,j,k=1 w l ijk |(U 1 ) ijk | + β 1 2 U 1 -X t -η t 1 2 F    (4.3b) X t+1 = argmin X β 0 2 U t+1 0 -A * X -η t 0 2 F + β 1 2 U t+1 1 -X -η t 1 2 F (4.3c) η t+1 0 = η t 0 -U t+1 0 + A * X t+1 (4.3d) η t+1 1 = η t 1 -U t+1 1 + X t+1 (4.3e)
The U 1 -subproblem (4.3b) is solved by soft-thresholding under nonnegative constraint. So the closed-form solution is given by

U t+1 1 ijk = max X t + η t 1 -w l ijk /β 1 , 0 .
The X -subproblem (4.3c) is a least squares problem. We rewrite the convolution into the componentwise multiplication by using the Fourier transform and the subproblem becomes

β 0 2 F{A} • F{X } -F U t+1 0 -η t 0 2 F + β 1 2 F{X } -F {U t+1 1 -η t 1 } 2 F .
Its closed-form solution reads as (4.4)

X t+1 = F -1 |F{A| 2 } + β 1 β 0 -1 F{A} • F{U t+1 0 -η t 0 } + β 1 β 0 F{U t+1 1 -η t 1 } ,
where |X | 2 and X are the componentwise operations of the square of absolute value of X , and complex conjugate of X , respectively. By assuming the boundary condition to be periodic, we use the 3D fast Fourier transformation (3D FFT) to compute (4.4) efficiently.

For the solution of U 0 -subproblem (4.3a), we need Theorem 4.1. The range of the following theorem goes beyond the solution of (4.3a). It gives a closed-form solution to a highly nonlinear and useful functional, namely the Thikhonov regularized KL divergence involving high-dimensional linear operators. It can also be seen as the closed-form solution of the proximal operator of this general KL divergence.

Theorem 4.1. Given β, b ∈ R 1 , η ∈ R m×n×d and G ∈ R m×n , consider the minimization problem (4.5) min U 1, T U -G log(T U + b 1) + β 2 U -η 2 F .
It has a closed-form solution

U = -B + B 2 + 4βC 2β ,
where

B = T T (1 + βb 1) -βη and C = T T (G -b 1) + βbT T T η.
Here T T is an operator from R m×n → R m×n×d with entries

(T T (K)) ijk = K ij , if k = d, 0, otherwise.
Proof. For any X ∈ R m×n×d and K ∈ R m×n , using the fact that T extracts the last slice, we easily deduce that

T (X ), K = m,n i,j=1 (T (X )) ij K ij = m,n i,j=1 X ijd K ij = X , T T K .
Therefore, T T defined as above is the adjoint operator for T . We directly get the derivative of 1, T U as

(4.6) ∇ 1, T U = ∇ T T 1, U = T T 1.
In addition, 1, G log(T U + b 1) can be rewritten as

m,n i,j=1 G ij log(U ijd + b 1). It is easy to get (4.7) ∂ ∂U uvw 1, G log(T U + b 1) = 0, if w = d.
When w = d, the partial derivative of 1, G log(T U + b 1) with respect to U uvw becomes

(4.8) ∂ ∂U uvw 1, G log(T U + b 1) = G uv U uvd + b .
Combining (4.7) and (4.8), we get

(4.9) ∇ 1, G log(T U + b) = T T   m,n i,j=1 G ij e T i T Ue j + b e i e T j   ,
where e i is the i-th canonical basis unit vector.

Denote

J (U) := 1, T U -G log(T U + b 1) + β 2 U -η 2 F . Combining (4.6
) and (4.9) as well as

∇ β 2 U -η 2 F = β(U -η), so ∇J (U) = T T 1 -T T   m,n i,j=1 G ij e T i T Ue j + b e i e T j   + β(U -η).
The minimizer of J with respect to U satisfies ∇J (U * ) = 0. We notice that U ijk with k = d are not involved in the sum between the parentheses and that the corresponding indexes in T are null, hence (4.10)

U * ijk = η ijk if k = d,
For those entries in the last slice of U * , we have

(4.11) 1 - G ij U * ijd + b + β(U * ijd -η ijk ) = 0, for all i = 1, • • • , m; j = 1, • • • , n
Multiplying by (U * ijd + b) both sides of (4.11), we have a quadratic equation with respect to each U ijd

β U * ijd 2 + (1 + βb -βη ijk )U * ijd -(G ij -b + βbη ijk ) = 0.
which is easily solved to yield (4.12)

U * ijd = -(1 + βb -βη ijk ) + (1 + βb -βη ijk ) 2 + 4β(G ij -b + βbη ijk ) 2β .
Combining (4.10) and (4.12) proves Theorem 4.1.

Based on Theorem 4.1, we have the following closed-form solution for U 0 -subproblem

U t+1 0 = -B + B 2 + 4β 0 C 2β 0 , where B = T T (1 + β 0 b) -µ(A * X t + η t 0 ) and C = T T (G -b 1) + β 0 bT T T (A * X t + η t 0 ).

4.2.

Removing false positives using a centroid method. For real data, point sources may be not on the grid, which means the discrete model is not accurate. In order to avoid missing sources, the regularization parameter µ is kept small, which leads to the over-fitting effect. That is to say that there are more point sources in our estimated result than the ground truth, which we must attempt to mitigate. Our optimization solution generally contains tightly clustered point sources, so we need to regard any such cluster of point sources as a single point source. The same phenomenon has been observed in [START_REF] Shuang | Generalized recovery algorithm for 3D super-resolution microscopy using rotating point spread functions[END_REF]. To this end, we apply a post-processing approach following [START_REF] Shuang | Generalized recovery algorithm for 3D super-resolution microscopy using rotating point spread functions[END_REF]. The method is based on the well-defined tolerance distance for recognizing clustered neighbors denoted as C. The criterion is to find nonzero entries that are within a certain distance from the given point whose pixel value is higher than or equal to other nonzero points in the neighborhood. For recognizing the cluster of point sources, we need to start from the entry with the highest intensity. By the above criterion, one clustered neighborhood is recognized. Then we compute the centroid of each cluster. That is to say, to get the centroid location (x, y, z) for one cluster C as

x = (i,j,k)∈C iX ijk (i,j,k)∈C X ijk ; y = (i,j,k)∈C jX ijk (i,j,k)∈C X ijk ; z = (i,j,k)∈C kX ijk (i,j,k)∈C X ijk .
The flux for this representative point source in cluster C is (i,j,k)∈C X ijk . We set the value of these entries in the recognized neighborhood as 0 and then we proceed with the searching process recursively. An estimated point source that cannot find out its corresponding ground truth is called false positive. To some false positives which may be due to the effect of the periodic boundary condition, we set a threshold, say 5% of highest intensity [START_REF] Shuang | Generalized recovery algorithm for 3D super-resolution microscopy using rotating point spread functions[END_REF]. Those entries whose pixel value is lower than the threshold will be regarded as false positives so that we directly delete them in the searching process. The whole post-processing method is summarized in Algorithm 4.2.

Algorithm 4.2 Removing cluster point sources by computing centroids (Centroid FP).

1: Set U as a zero 3D tensor of the same size as X ; 2: Find the maximum value of the solution X ; 3: Recognize the neighbour C based on the above criterion and recursively check the other nonzero points; 4: Compute the centroid of C and set the value of this centroid entry in U as the summation of all the pixel values in C; 5: Set X (C) = 0; 6: If there is any nonzero entry in X , go to Step 2, otherwise, go to the next step; 7: Set pixel value of these entries in U whose value is lower than 5% of the highest intensity as 0.

Estimating the flux values.

In Section 4, we provided an algorithm for estimating the locations and fluxes for the point sources. However, numerical results show that the flux values are generally underestimated. In [START_REF] Shuang | Generalized recovery algorithm for 3D super-resolution microscopy using rotating point spread functions[END_REF][START_REF] Min | FALCON: fast and unbiased reconstruction of high-density super-resolution microscopy data[END_REF], least squares fitting is used for improving the resolution as well as updating the corresponding fluxes. However, our problem is not Gaussian noise or additive Poisson noise as used in their paper. Our Poisson noise is datadependent, which cannot force the regenerated image T (A * X * ) to match the observed data with least squares, where X * is the result after solving the nonconvex optimization problem and post-processing. Our aim is thus to estimate the source fluxes from the KL data-fitting term appropriate to the Poisson noise model when the source 3D positions have already been accurately estimated.

Let the PSF corresponding to the i-th source be arranged as the column vector h i . The stacking of the P column vectors in the same sequence as the source labels for the P sources then defines a system PSF matrix H, with H = [h 1 , h 2 , • • • , h P ] ∈ R K×P , where K is the total number of pixels in the vectorized data array, so K = mn. The vectorized observed image is denoted by g ∈ R K×1 . The uniform background is denoted as the vector b1 with 1 ∈ R K×1 . The flux vector is denoted as f ∈ R P ×1 . Here the problem is overdetermined meaning that the number of point sources P is much smaller than the number of available data K = mn. Therefore we need to do some refinement of the estimates by minimizing directly data fitting term. Since the negative log-likelihood function for the Poisson model, up to certain data dependent terms, is simply the KL divergence function, D KL (Hf + b1, g) = 1, Hf -g log(Hf + b1) , its minimization with respect to the flux vector f , performed by setting the gradient of D KL with respect to f (see [START_REF] Harmany | This is spiral-tap: Sparse poisson intensity reconstruction algorithmstheory and practice[END_REF]) zero, yields the nonlinear relation

∇D KL (Hf + b1, g) =H T 1 - m i=1 g i e T i (Hf + b1) H T e i = m i=1
e T i (Hf + b1 -g) e T i (Hf + b1)

H T e i = 0.

(4.13)

Consider now an iterative solution of (4.13). Since b is proportional to the vector of ones, 1, we write (4.13) as 0 = bD KL (Hf + b1, g)

= K i=1 e T i (Hf + b1 -g) e T i (b1) e T i (Hf + b1) H T e i = K i=1 e T i (Hf + b1 -g) e T i (Hf + b1 -Hf ) e T i (Hf + b1) H T e i = K i=1 e T i (Hf + b1 -g) H T e i - K i=1 e T i (Hf + b1 -g) e T i Hf e T i (Hf + b1) H T e i = H T Hf + H T (b1 -g) - K i=1 e T i (Hf + b1 -g) e T
i Hf e T i (Hf + b1)

H T e i . In this form, we easily get

H T Hf = H T (g -b1) + K i=1 e T i (Hf + b1 -g) e T i Hf e T i (Hf + b1) H T e i ,
By multiplying by (H T H) -1 the two sides of this equation, we directly get

(4.14) f = f G + K i=1 e T i (Hf + b1 -g) e T
i Hf e T i (Hf + b1)

H + e i ,
where

H + = (H T H) -1 H T and f G = H + (g -b1)
is the solution corresponding to the Gaussian noise model. This suggests the following iterative algorithm:

(4.15)

f n+1 = f G + K(f n ), n = 1, 2, • • •
where

K(f ) = K i=1 e T i (Hf + b 1 -g) e T
i Hf e T i (Hf + b1)

H + e i .

To emphasize that our nonconvex optimization model is based on the use of a KL data fitting (KL) term and a nonconvex (NC) regularization term, we designate our approach as KL-NC that has been enhanced with post-processing and a refined estimation of flux. We summarize the different steps of our proposed method for the rotating PSF problem in Algorithm 4.3.

Numerical experiments.

In this section we apply our optimization approach to solving simulated rotating PSF problems for point source localization and compare it to some other optimization methods. The codes of our algorithm and the others with which we compared our method were written in MATLAB 9.0 (R2016a), and all the numerical experiments were conducted on a typical personal computer with a standard CPU (Intel i7-6700, 3.4GHz).

The fidelity of localization is assessed in terms of the recall rate, defined as the ratio of the number of identified true positive point sources over the number of true positive point sources, and the precision rate, defined as the ratio of the number of identified true positive point sources over the number of all point sources obtained by the algorithm; see [START_REF] Birk | Super-Resolution Microscopy: A Practical Guide[END_REF].

To distinguish true positives from false positives from the estimated point sources, we need to determine the minimum total distance between the estimated point sources and true point sources. Here all 2D simulated observed images are described by 96-by-96 matrices. We set the number of zones of the spiral phase mask responsible for the rotating PSF at L = 7 and the aperture-plane side length as 4 which sets the pixel resolution in the 2D image (FFT) plane as 1/4 in units of λz I /R. The dictionary corresponding to our discretized 3D space contains 21 slices in the axial direction, with the corresponding values of the defocus parameter, ζ, distributed uniformly over the range, [-21, 21]. According to the Abbe-Rayleigh resolution criterion, two point sources that are within (1/2)λz I /R of each other and lying in the same transverse plane cannot be separated in the limit of low intensities. In view of this criterion and our choice of the aperture-plane side length and if we assume conservatively that our algorithm does not yield any significant superresolution, we must regard two point sources that are within 2 image pixel units of each other as a single point source. Analogously, two sources along the same line of sight (i.e., with the same x, y coordinates) that are axially separated from each other within a single unit of ζ must also be regarded as a single point source.

As for real problems, our simulation does not assume that the point sources are on the grid points. Rather, a number of point sources are randomly generated in a 3D continuous image space with certain fluxes. We consider a variety of source densities, from 5 point sources to 40 point sources in the same size space. For each density, we randomly generate 20 observed images and use them for training the parameters in our algorithm, and then test 50 simulated images with the well-selected parameters. The number of photons emitted by each point source follows a Poisson distribution with mean of 2000 photons. Instead of adding Poisson noise as additive noise as in [START_REF] Shuang | Generalized recovery algorithm for 3D super-resolution microscopy using rotating point spread functions[END_REF], we apply data-dependent Poisson noise by using the MATLAB command

G = poissrnd(I0+b),
where I0 is the 2D original image formed by adding all the images of the point sources, and b is the background noise which we set to a typical value 5. Here, poissrnd is the MATLAB command whose input is the mean of the Poisson distribution.

3D localizations for low and high density cases.

In this subsection, we will show some 3D localizations results in low and high density cases (15 point sources and 30 point sources); see Figure 3.

Note that our 3D localization is estimated very well for 15 point sources. In this example, all the true positives are identified and there is only one false positive. From Figure 3(c) and (e), we see there are two true point sources that are close in their transverse coordinates, but well-separated in axial direction. In the 2D original image (without noise) in Figure 3(c), the corresponding two rotating PSF images are not overlapping and our approach can well identify both two point sources.

Figure 3(b), (d) and (f) show the high density (30 point sources) case with many overlapping rotating PSF images. Such image overlap in the presence of Poisson noise makes the problem difficult. The number of point sources is not easily obtained by observation. In this specific case, our algorithm still identifies all the true point sources correctly, but produces 9 false positives. From Figure 3(d), we can see that these false positives come from the serious PSF overlapping. Recall that the solution of our nonconvex optimization problem involves outer and inner iterations. Here, we set the maximum number of outer iterations as 2 and maximum number of inner iterations as 400. From Figures 4 and5, we can see how and why we chose these two maximum iteration numbers. Let us denote the number of maximum outer iterations and inner iterations as Max out and Max in respectively. In Figure 4, we fixed Max in as 400 and tried different Max out . As we see, the recall rate, whether with or without post-processing, changes little. The precision rate increases, however, when we change Max out from 1 to 2, but it changes considerably less when changing Max out from 2 to 3, while the computational time increases linearly with Max out . In Figure 5, on the other hand, we fixed Max out as 2 and changed Max in . Here Max in = 400 gives the best results with a relatively low time cost. In both cases, the precision rate can be improved significantly with post-processing. Here, we consider the practical convergence of our algorithm. That is to check the differences between the auxiliary variables (U 0 and U 1 ) and their corresponding substituted variables (A * X and X ), i.e., the value of U 0 -A * X F and U 1 -X F in each iteration. Noted that the maximum number of inner and outer iterations are 400 and 2 respectively. We plot the values of these two terms for each iteration in Figures 6 and7. Both two terms decrease iteratively. Note that the values of U 0 F and A * X F are about 3000 and the values of U 1 F and X are about 400. Therefore, the differences between auxiliary variables and their substituted variables are small, which shows the convergence numerically.

Comparison of Algorithm KL-NC with other algorithms.

In this subsection we compare our approach with three other optimization algorithms: KL-1 (KL data fitting with In KL-1 , we need to solve the following optimization problem min

X ≥0 {D KL (T (A * X ) + b 1, G) + µ X 1 } .
We use ADMM by introducing two auxiliary variables U 0 = A * X and U 1 = X . We get a similar augmented Lagrangrian function with w l ijk replaced by µ in (4.2). Therefore, we use Theorem 4.1 to solve this problem. In 2 -1 , we need to solve the following optimization problem min

X ≥0 1 2 T (A * X ) + b 1 -G 2 F + µ X 1 .
The problem is similar to the optimization problem in [START_REF] Shuang | Generalized recovery algorithm for 3D super-resolution microscopy using rotating point spread functions[END_REF], except that we include the background which is estimated easily by signal processing. We can therefore still use ADMM for solving this problem. In 2 -NC, we need to solve the following optimization problem min

X ≥0 1 2 T (A * X ) + b 1 -G 2 F + µR(X ) ,
where R(X ) is the same as KL-NC method. We see that IRL1 can also work well for solving this problem with the 1 weighted optimization as a subproblem which we solve by ADMM since it is similar to 2 -1 .

Both the initial guesses of X and U 0 are set as 0 for all these methods. In order to do the comparison, we plot the solutions for the four models as well as the ground truth in the same space; see Figures 8 and9 which correspond to the cases of 15 and 30 point sources, respectively.

From Figure 8 and Figure 9, we see the overfitting of the 1 regularization models (KL-1 and 2 -1 ). Before post-processing, the solutions of these two algorithms spread out the PSFs a lot and have many false positives. After post-processing, both algorithms are improved, especially KL-1 . However, in comparison to the nonconvex regularization (KL-NC and 1 -NC), they still have many more false positives. Among the four algorithms, our approach (KL-NC) performs the best in terms of the recall and precision rates.

We also tested a number of other values of point source density, namely 5, 10, 15, 20, 30 and 40 point sources, and computed the average recall and precision rate of 50 images for each density and for each algorithm; see Tables 1 and2. The results show the superior results of our method in terms of both recall and precision rates. In Table 1, the best recall and precision rates in each case are labeled by bold fonts. As in the above discussion, our nonconvex regularization tends to eliminate more false positives which increases the precision rate. The KL data-fitting term, on the other hand, improves the recall rate as we see by comparing the results of KL-NC with 2 -NC. Before post-processing, we see that all the algorithms have low precision rates, especially the two employing the 1 regularization model at less than 10%.

In the following, we compare the results of the estimation of the flux in these four algorithms, considering specifically the case of 15 point sources. In Figure 10, we plot the fluxes of ground truth as well as the fluxes of the estimated point sources for the true positive point sources. For the false positive point sources, we only show the estimated fluxes. Both 1 models underestimate the fluxes. The rotating PSF images for false positives carry the energy away from the true positive source fluxes. In nonconvex models, we also have similar observations when we have false positives. For example, in Figure 10(b), we see the flux on the third bar is underestimated more than the others. We note that it's rotating PSF is overlapping with the image of a false positive. The more false positives an algorithm recovers the more they will spread out the intensity, leading to more underestimated fluxes for the true positives. We also tested 50 different observed images for each density, and analyzed the relative where (f est , f tru ) is the pair which contains the flux of an identified true positive and the corresponsing ground truth flux. We plot the histogram of the relative errors on these four models in Figures 11 and12. We still see the advantage of KL-NC over other algorithms in this respect. The distribution of relative errors mostly lies within [0, 0.1]. For the 1 regularization algorithms, the distribution of the relative error spreads out and there are many cases with error higher than 0.3. For 2 -NC, we have the similar result for the 15 point sources case; however, in higher density cases, KL-NC has better results comparing with 2 -NC.

Low number of photons case.

In the previous subsection, we saw the advantages of using our nonconconvex regularization term. However, KL-NC does not seem to possess significant advantages over the 2 -NC model, especially in the cases of low source density. This is because only when we have a sufficiently large number of photons, as was the case in the results discussed earlier, the 2 data-fitting term corresponding to the maximum likelihood estimator for the Gaussian noise model is well approximated by the KL data-fitting term for the Poisson noise model assumed here. In order to verify this, we study the case of low photon number in this subsection. In Figure 13, we see the observed image with three different photon numbers per source, namely 500, 1000, and 2000. For the case of 500 photons per source, the image is very dim, and we can barely distinguish the rotating PSF signal from the background noise, so we do not consider this case any further. For the case of 1000 photons being emitted by each point source, the corresponding image follows a Poisson distribution with a mean of 1000 photons. We take each image to contain 15 point sources. As before, we first randomly generated 20 images, which were used for training parameters. Then we tested 50 different observed images by using the trained parameter values. We got 90.53% for the recall rate and 78.32% for the precision rate in the KL-NC algorithm, but only 81.33% and 64.24%, respectively, for the two rates in the 2 -NC algorithm. For evaluation of the flux estimates, we refer to Figure 14 from which we see that once again KL-NC gives much better results than 2 -NC since the bars corresponding to the low relative error case are higher. We therefore conclude that the KL data-fitting term is much better than the 2 data-fitting term when dealing with low numbers of photons. 6. Conclusions and Current Directions. In this paper, we have proposed a nonconvex optimization algorithm for the 3D localization of a swarm of randomly spaced point sources using a specific rotating PSF which has a single lobe in the image of each point source. Rotating PSFs that have a single lobe have several obvious advantages over those that have multiple lobes or are more spread out, particularly when dealing with relatively high source densities and relatively small photon numbers per source. We employed a post-processing step based both on centroiding the locations of recovered sources that are tightly clustered and thresholding the recovered flux values to eliminate obvious false positives from our recovery sources. In addition, we have proposed a new iterative scheme for refining the estimate of the source fluxes after the sources have been localized. These techniques can be applied to other rotating PSFs as well as other depth-encoding PSFs for accurate 3D localization and flux recovery of point sources in a scene from its image data under the Poisson noise model. Applications include not only 3D localization of space debris, but also super-resolution 3D single-molecule localization microscopy, e.g. [START_REF] Birk | Super-Resolution Microscopy: A Practical Guide[END_REF][START_REF] Von Diezmann | Three-dimensional localization of single molecules for super-resolution imaging and single-particle tracking[END_REF]. Tests of this algorithm based on real data collected using phase masks fabricated for both applications are currently being planned. In addition, work involving snapshot multi-spectral imaging, which will permit accurate material characterization, as well as higher 3D resolution and localization of space microdebris via a sequence of snapshots is under way.

Figure 1 .

 1 Figure 1. Illustration of the discretized forward model for the single lobe PSF.

Figure 2 .

 2 Figure 2. Plot of θ(t; a) and θ (t; a). In (a), we choose different values of a; in (b) we choose different values of while fixing a = 80.

Algorithm 4 . 3 1 : 2 : 3 :

 43123 KL-NC with post-processing and estimation of the flux. Input: X 0 ∈ R m×n×d and G ∈ R m×n . Output: The locations of point sources with the corresponding flux {(x * i , y * i , z * i , f * i )}. Solve the nonconvex optimization problem by Algorithm 4.1 and get the minimizer X * ; Do the post-processing by Algorithm 4.2 to get the locations of estimated point sources {(x * i , y * i , z * i )}; Estimate the fluxes by the iterative scheme (4.15) and get f * .

Figure 3 .

 3 Figure 3. 3D localizations for low and high density cases: (a), (c) and (e) are the results for the 15 point sources case; (b), (d) and (f ) are the results for the 30 point sources case. "o" is the ground truth point source and "+" is the estimated point source.

Figure 4 .

 4 Figure 4. Tests on effects of the maximum number of outer iterations.

Figure 5 .Figure 6 .

 56 Figure 5. Tests on effects of the maximum number of inner iterations.

Figure 7 .

 7 Figure 7. U0 -A * X F in the each iterations.

Figure 8 .

 8 Figure 8. Solution X from 4 algorithms (15 point sources).

Figure 9 .

 9 Figure 9. Solution X from 4 algorithms (30 point sources).

Figure 10 .

 10 Figure 10. Tests on estimating flux values. The bar graph with no value ground truth part corresponds to a false positive.

Figure 11 .

 11 Figure 11. Histogram of the relative errors of the flux values for the 15 point sources case.

Figure 12 .

 12 Figure 12. Histogram of relative errors of flux values in the 30 point sources case.

Figure 13 .Figure 14 .

 1314 Figure 13. Observed image with different numbers of photons in each point source.

Table 1

 1 Comparison of KL-NC with KL-1, 2-1 and 2-NC. All the results are after performing postprocessing.

		KL-1	KL-NC	2 -1		2 -NC
	No.source Recall	Prec.	Recall	Prec.	Recall	Prec.	Recall	Prec.
	5	98.53% 58.67% 100.00% 93.71% 100.00% 68.41% 97.20% 87.95%
	10	99.40% 65.33% 100.00% 92.83%	99.60%	56.54% 94.40% 83.33%
	15	98.53% 58.67%	98.80%	87.22% 99.20% 54.20% 93.33% 83.45%
	20	98.10% 57.76% 98.60% 85.29%	97.70%	56.03% 93.60% 78.89%
	30	95.33% 54.54%	94.47%	79.93% 95.67% 54.16% 89.73% 74.63%
	40	93.60% 53.87% 95.40% 69.22%	95.25%	49.32% 82.45% 51.63%

Table 2

 2 Comparison of KL-NC with KL-1, 2-1 and 2-NC. All the results are without post-processing.

		KL-1	KL-NC	2 -1		2 -NC
	No.source	Recall	Prec.	Recall	Prec.	Recall	Prec.	Recall	Prec.
	5	100.00% 2.82% 100.00% 44.85% 100.00% 8.89% 98.80% 57.67%
	10	100.00% 2.91% 100.00% 44.87% 100.00% 8.15% 95.60% 55.45%
	15	100.00% 2.82% 100.00% 44.67% 100.00% 5.95% 96.80% 56.74%
	20	100.00% 3.26% 99.70% 46.29% 100.00% 4.83% 96.20% 52.15%
	30	100.00% 3.56% 98.27% 44.44% 100.00% 4.85% 95.07% 52.02%
	40	100.00% 3.75% 99.00% 23.42% 100.00% 6.07% 98.25% 32.19%
	error in the estimated flux values, which we define as
				error =	f est -f tru f tru	,

regularization model);

-1 (least squares fitting term with 1 regularization model) and 2 -NC (least squares fitting term with nonconvex regularization model). For all these models, we do the same post-processing and estimation of flux values after solving the corresponding optimization problem.

The work of the second author was supported by HKRGC Grants No. CUHK14306316, HKRGC CRF Grant C1007-15G, HKRGC AoE Grant AoE/M-05/12, CUHK DAG No. 4053211, and CUHK FIS Grant No. 1907303. The work of the third author was partially funded by the French Research Agency (ANR) under grant No ANR-14-CE27-001 (MIRIAM) and by the Isaac Newton Institute for Mathematical Sciences for support and hospitality during the programme Variational Methods and Effective Algorithms for Imaging and Vision, EPSRC grant no EP/K032208/1. The work of the fourth author was also supported by HKRGC Grants No. CUHK14306316. The research of the first, fourth and fifth author were supported by the US Air Force Office of Scientific Research under grant FA9550-15-1-0286.