A content-adaptive unstructured grid based integral equation method with the TV regularization for SPECT reconstruction

Existing reconstruction methods for single photon emission computed tomography (SPECT) are most based on discrete models, leading to low accuracy in reconstruction. Reconstruction methods based on integral equation models (IEMs) with a higher order piecewise polynomial discretization on the pixel grid for SEPCT imaging were recently proposed to overcome the accuracy deficiency of the discrete models. Discretization of IEMs based on the pixel grid leads to a system of a large dimension, which may require higher computational costs to solve. We develop a SPECT reconstruction method which employs an IEM of the SPECT data acquisition process and discretizes it on a content-adaptive unstructured grid (CAUG) with the total variation (TV) regularization aiming at reducing computational costs of the integral equation method. Specifically, we design a CAUG of the image domain for the discretization of the IEM, and propose a TV regularization defined on the CAUG for the resulting ill-posed problem. We then apply a preconditioned fixed-point proximity algorithm to solve the resulting non-smooth optimization problem, and provide convergence analysis of the algorithm. Numerical experiments are presented to demonstrate the superiority of the proposed method over the competing methods in terms of suppressing noise, preserving edges and reducing computational costs.

method. Specifically, we design a CAUG of the image domain for the discretization of the IEM, and propose a TV regularization defined on the CAUG for the resulting ill-posed problem. We then apply a preconditioned fixed-point proximity algorithm to solve the resulting non-smooth optimization problem, and provide convergence analysis of the algorithm. Numerical experiments are presented to demonstrate the superiority of the proposed method over the competing methods in terms of suppressing noise, preserving edges and reducing computational costs.
1. Introduction. Single photon emission computed tomography (SPECT) can visually provide patients' functional information by estimating distribution of radioactive tracer inside their bodies through reconstruction of measured projection data [30]. Mathematically, SPECT reconstruction requires solving an ill-posed problem. Conventional reconstruction methods widely use discrete models, piecewise constant approximations of the integral equation model describing the SPECT data acquisition process, which lead to low accuracy in reconstruction. To overcome the accuracy deficiency of the conventional reconstruction methods based on the discrete models, reconstruction methods based on the integral equation model with discretization using higher order piecewise polynomials defined on the pixel grid were recently proposed in [11]. The use of the pixel grid leads to discrete systems of large dimensions and requires more computing time to solve them. To address this issue, we propose a content-adaptive unstructured grid (CAUG) based integral equation method with the total variation (TV) regularization for SPECT reconstruction.
SPECT reconstruction has been a long standing research topic [3,10,11,12,13,19,20,23,28]. Existing methods were most based on SPECT discrete models due to their convenience in implementation and consistence with usual sampling methods [18]. These models are in fact certain piecewise constant approximations of the continuous model in an integral equation form, which describes the SPECT data acquisition process. Due to its low accuracy, the piecewise constant approximation introduces a model error that can not be compensated by any techniques employed in processing the resulting discrete system for the reconstruction [16,18]. To reduce the model error, recent paper [11] developed an integral equation method by employing a continuous model considering the SPECT system spatial resolution and photon attenuation, constructing pixel-based piecewise polynomial approximation, regularizing the resulting discrete system by penalizing discontinuities and noise of the reconstructed images. A new challenge of the reconstruction method developed in [11] lies in its large dimension of the resulting discrete system due to the use of the pixel grid in the discretization of the integral equation model. The large dimension of such a system transfers to high computational costs in solving it.
In passing, we point out that reconstruction methods based on continuous models with discretization on a CAUG were previously studied in the literature [2,3] for tomographic reconstruction. However, the absent of appropriate regularization in these methods prevents them from effectively overcoming the ill-posedness of the resulting discrete systems. The method proposed in [2] did not use regularization while in [3] a Gibbs prior using the potential function was applied to tomographic reconstruction on the CAUG. Using such a function as a regularization fails to effectively suppress noise with preserving important features of the reconstructed images due to smoothness of the function. As a result, the reconstructed images produced by these methods suffer from low quality due to their contamination of noise or/and the present of staircase artifacts. This paper addresses the large dimensionality challenge of the integral equation method proposed in [11] by employing a CAUG for the discretization of the SPECT integral equation model aiming at reducing its computational complexity while preserving its nice features such as high accuracy and computational stability. Specially, we consider an integral equation which models the SPECT photon attenuation. As a step of proving the concept, we consider a simple version of the model that does not address system spatial resolution issue. A more comprehensive model will be considered in a future occasion. We design a CAUG for the image domain based on an initial estimation of the underlying image and discretize the integral equation through a collocation method with piecewise linear basis functions defined on the CAUG. We further develop the TV regularization method defined on the CAUG for the resulting ill-posed linear system, aiming at suppressing noise and preserving important features of the reconstructed images. The regularized SPECT reconstruction problem is then formulated as a non-smooth convex optimization problem via the penalized maximum likelihood criterion. A preconditioned fixed-point proximity algorithm is applied to address the non-smoothness of the optimization problem and to solve it. The algorithm inherits the advantages of the fixed-point proximity algorithm and the maximum-likelihood expectationmaximization (ML-EM) algorithm. We provide rigourous analysis for convergence of the proposed algorithm. Numerical experiments are presented to demonstrate the accuracy and effectiveness of the developed method for SPECT reconstruction. We compare the proposed method with the ML-EM and the maximum a posteriori EM (MAP-EM) on the CAUG, and the traditional discrete TV regularized reconstruction method in the three aspects of suppressing noise, preserving edges and reducing computational costs. Numerical results show that the developed method is superior to the competing methods in all these accounts. This paper is organized in nine sections. In section 2, we describe the integral equation model for SPECT to be used in this paper. We design in section 3 a CAUG for the image domain and introduce basis functions based on the CAUG for the representation of the underlying image. We apply in section 4 the collocation method using the basis functions constructed in section 3 to discretize the integral equation. Section 5 is devoted to presenting an edge-preserving regularization method defined on the CAUG for the resulting ill-posed problem, and formulating a non-smooth convex optimization problem via the penalized ML criterion. In section 6, we employ a preconditioned fixed-point proximity algorithm to solve the resulting optimization problem. We provide in section 7 a rigourous convergence analysis for the iterative algorithm. In section 8 we present numerical experiments to compare the proposed method with several competing methods for SPECT reconstruction. We conclude this paper in section 9.
2. The integral equation model for SPECT. We describe in this section an integral equation for SPECT which models the photon attenuation of the SPECT imaging system.
SPECT seeks a function (image) of the radioactivity distribution of a phantom from given projection data. The underlying projection data acquisition process may be naturally formulated as an integral equation model, involving a kernel describing photon attenuation which obeys Beer's law. In particular, let g : Y × Θ → R be the continuous distribution function on the projection domain with Y ⊂ R as an interval on the detector face and Θ := [0, 2π] as the angular rotation range, and f be the radiotracer distribution on the square image domain Λ ⊂ R 2 . Given an attenuation distribution µ : Λ → R, we define the kernel as and further define the integral operator K : C(Λ) → C(Y ×Θ) in terms of the kernel K as Here L(ν, θ) denotes the path of the collimated gamma ray corresponding to the lateral position ν on the detector at projection angle θ, the line segment L(ν, θ, x) is the attenuation path along L(ν, θ) from the spatial position x := (x, y) ∈ Λ to the detector, dl is the arc length differential, and C(Υ) is the space of all continuous functions on Υ. As a result, we calculate the projection data g of f and µ along a single line by the following first-kind Fredholm integral equation Note that the operator K is the well-studied attenuated Radon transform [30]. Clearly, the integral equation (1) is an ill-posed problem due to smoothness of the kernel, and hence, its solution does not continuously depend on the given projection data.
Conventional methods for solving the integral equation model (1) represent the underlying function and the integral operator directly via pixel-based piecewise constant basis functions. This leads to a discrete model for SPECT that suffers from the model error and heavy computational burden. To overcome these difficulties, we shall propose a continuous representation of the underlying function based on a CAUG, and solve the integral equation model using the collocation method with continuous basis functions defined on the CAUG. The use of the CAUG provides a sparse and compact image representation, and substantially reduces the number of spatial samples on the image domain. This in turn leads to fast numerical computation.
3. Content-adaptive unstructured grid. We design in this section a CAUG for the image domain and describe the related basis functions on the CAUG. We first recall the description of an unstructured grid for the image domain. We then construct a CAUG for the image domain using the quadtree scheme and the force equilibrium method. We further introduce basis functions on the CAUG.
We first describe an unstructured grid for the image domain. Let N c := {1, 2, . . . , c} be the set of positive integers for a positive integer c in this paper. An unstructured grid C := {C 1 , C 2 , . . . , C q } of the image domain Λ is a set of q triangles, having the following properties: vertex or a common edge of C i and C j . (iii) Each triangle can be indicated by listing their vertices, in the anti-clockwise direction.
We next design a CAUG for the image domain Λ from a given estimate f 0 (e.g., figure 1(a)) of the underlying function. The estimate may be reconstructed by a pixel-based efficient iterative algorithm from measured projection data. We generate the CAUG for the image domain by the estimate and the given attenuation distribution through the following two steps.
The first step yields the quadtree grid G 0 by the quadtree scheme [26] to obtain a set of vertices. For adaptive algorithms of this context, see [6,33]. For the description of the quadtree scheme, we need a local error function defined for a function h 0 : Λ → R and a triangle ∆ ∈ G 0 by Here, Π 1 denotes the space of polynomials of degree at most 1 on R 2 . With a threshold 1 > 0, the quadtree scheme can be described briefly as follows: (i) Input a function h 0 : Λ → R, and construct an initial quadtree grid G 0 = {∆ 1 , ∆ 2 } through splitting the image domain Λ into two triangles ∆ 1 and ∆ 2 via one of its diagonals.
if E(h 0 , ∆ 0 ) < 1 , stop and output the quadtree grid G 0 (h 0 ). (iii) Otherwise, update the grid G 0 by splitting ∆ 0 into four triangles self-similar to itself. Goto (ii). We first apply the above quadtree scheme to the initial estimate f 0 and obtain the grid G 0 (f 0 ), which is adaptive to the structure of the underlying function carried by f 0 . Since the attenuation distribution provides anatomical information of a phantom, we may improve the grid G 0 (f 0 ) through employing the quadtree scheme to the given attenuation distribution µ (e.g., figure 1(b)) with G 0 (f 0 ) as the initial quadtree grid. This process results in a quadtree grid G 0 (µ) adapted to both f 0 and µ. By the estimate f 0 and the distribution µ, we may divide Λ into two parts: the content region Ω and the background region Λ \ Ω, where the content region is the region with radioactive tracer and the attenuation distribution (e.g., the large elliptic region in figures 1(a)-(b)). We remove the triangles and the related vertices outside the content region. In this way, we obtain the resulting quadtree grid G 0 (figure 1(c)) of Ω by restricting the quadtree grid G 0 (µ) of Λ to the content region Ω. Clearly, we have that Let V 0 := {v 0 j = (x j , y j ) : j ∈ N d } denote the set of vertices of G 0 . The second step applies the force equilibrium method [24] to the vertex set V 0 yielded by f 0 to generate the CAUG of the domain Ω. We seek the CAUG whose each element is "equilateral" for a more accurate discretization of (1). As indicated in [24], each element of the grid generated by the force equilibrium method tends to be equilateral. The method consists of the Delaunay triangulation [8] of the vertex set and its equilibration. Using the former we can generate an initial CAUG whose elements may not be equilateral, and using the latter we may yield a set of equilibrium vertices from the initial grid. The resulting set leads to the CAUG to be sought. To get such a set, the method considers a simple mechanical analogy between an unstructured grid and a structure of springs. In this model, each edge of the grid corresponds to a spring, and each vertex corresponds to a joint. Each edge has a force-displacement relationship depending on its current length and its desired length, leading to the movement of the connected vertices. Moreover, in this paper we set a weight for each vertex based on density of the vertex set to prevent the vertex from moving. Using the relationship and the weight that is determined by f 0 , we obtain a desired CAUG which can efficiently fit the given estimate f 0 .
We describe the force equilibrium method based on the above introduction. For an initial CAUG yielded by the Delaunay triangulation from V 0 , we assume that it has n 1 edges. Using the l-th (l ∈ N n1 ) edge with its current length l and its two endpoints v 0 i , v 0 j ∈ V 0 , by [24] we may define the function characterizing the force-displacement relationship acting on its endpoints as denote the function that models the elastic force acting on v 0 j , where R j denotes the index set of edges connected to v 0 j . We further let m j be the weight for v 0 j which is proportional to 1/S j , where S j is the mean of areas of all triangles having v 0 j . Using the weight, we define the function that models the frictional force of v 0 j as We then have the following force function on v 0 j , given by Based on the above definitions, the force equilibrium method yields an unstructured grid whose elements tend to be equilateral. In terms of vertices, the method finds a set of equilibrium vertices by solving for the following system where V is an d×2 matrix. A simple approach to solve (6) is to introduce an artificial time-dependence [24]. At the discretized (artificial) time t k = k∆t (k = 0, 1, 2, . . .), the approximate solution V k is updated by where V 0 is the d×2 matrix whose first column contains x-coordinates of all vertices in V 0 and second column contains y-coordinates of these vertices. With a threshold 2 > 0, the force equilibrium method can be described briefly as follows: (i) Input the vertex set V 0 , obtain the matrix V 0 and yield an initial CAUG G 1 from the set through the Delaunay triangulation.
stop and output the grid.
(iii) Otherwise, the grid G 1 is iteratively improved by (7); update the vertices set V 0 . Goto (i). Due to the convexity of the yielded grid, we need to recover edges of Ω by removing triangles whose centroids are not in the content region. We then obtain the resulting CAUG (figure 1(d)) for the domain Ω which consists of n triangles, denoted by G := {∆ i : i ∈ N n }. We further have the set of vertices of G, defined by V : The quadtree grid G 0 generated by the quadtree scheme from f 0 and µ; (d) The CAUG from G 0 through the force equilibrium method.
We then introduce basis functions defined on the generated CAUG for the representation of the underlying image. Some definitions are required to describe that. For the resulting grid G, these triangles are called grid elements, indicated by listing their vertices in the anti-clockwise direction. This implies that for G and V, there exists a relationship between their index sets N n and N d , characterized as an n × 3 matrix N with N (i, k) ∈ N d (i ∈ N n , k ∈ N 3 ) as the index of the k-th vertex of the i-th grid element. Let v ik := v N (i,k) , we may define the set T of grid elements of G as (8) T : When using a continuous piecewise linear function for G, by the set (8) the linear function within each grid element is uniquely determined by its values at the vertices of the element [29]. Hence, by the vertex set V we may define piecewise linear basis functions as We obtain a finite-dimensional operator equation through applying the collocation method to solve the integral equation (1). The collocation method [5] provides an approach to approximate a solution of the integral equation by a sequence of functions in finite-dimensional subspaces of C(R 2 ). The method has lower computational cost and uses linear functionals. Let {l t,i : t ∈ N N1 , i ∈ N N2 } be the set of linear functionals, defined for a function ρ ∈ C(Y × Θ) by where i is the interval of the i-th detector bin, θ t denotes the t-th projection angle, N 1 is the number of projection angles and N 2 is the number of detector bin. Let X d be an d-dimensional subspace of C(R 2 ), spanned by the basis {e j : j ∈ N d }. With the above analyses, the collocation method for solving the integral equation (1) is to find a solutionf ∈ X d such that We now transform the finite-dimensional operator equation (9) to a linear system using the defined basis. Sincef ∈ X d , the solutionf is represented by the following formula (10) where f j :=f (x j , y j ) is the representation coefficient of the functionf at the j-th basis (i.e., the estimate off at the j-th vertex). Since the operator K is linear, by (10) we find that for the left-hand side of (9), where Ke j is the footprint function of the j-th basis [31]. For each t ∈ N N1 , we let where (A t ) ij denotes the value of projecting the j-th basis onto the i-th detector bin and (g t ) i is the expected projection data received by the i-th detector bin at projection angle θ t . With (11)-(12), we then obtain a compact formula ] is the expected projection data at θ t . Moreover, by accumulating (13) at each angle in column wise we let and find that the integral equation (1) can be approximated by the following linear system Here m := N 1 × N 2 , g ∈ R m is the expected projection data at all angles, and A ∈ R m×d is the overall system matrix. In a SPECT system, the photon detection actually follows a temporal Poisson distribution. We consider the vector r ∈ R m containing expected counts due to background activity. In consideration of background activity and the randomness in the projection data, the measured projection data g ∈ R m , which relate to the representation coefficients f through the expectation vector g and the linear system (15), can be described by the Poisson model (16) Poisson(Af + r) = g.
The goal of SPECT reconstruction is to estimate f from a given realization of g so that the resulting solutionf preserves image features as much as possible, with suppressing Poisson noise. To this end, regularization is necessary to solve (16) due to the ill-posedness of the problem. However, only the potential function as a regularization was used to tomographic reconstruction on the CAUG. Due to smoothness of the potential function, it may not be effective for improving the quality of the reconstructed images.

5.
Regularization defined on the CAUG. We study a regularization method and an optimization model based on the CAUG in this section. We first propose a regularization method defined on the CAUG from the TV regularization. We then formulate the regularized SPECT reconstruction problem as an optimization model via the penalized ML criterion for solving the Poisson model.
We propose a regularization method defined on the CAUG from the TV [25,27], aiming at suppressing noise and preserving important features of the reconstructed images. In the context of the pixel grid, regularization methods have been intensively studied, especially the piecewise constant approximation [22,23] of the anisotropic TV. This discretization of the anisotropic TV is effective to these categories due to the anisotropic structures of images. Moreover, it can be characterized as the composition of the 1 -norm and the first-order difference matrix. The 1 -norm has a proximity operator with a closed form, and the composite structure admits a fixed-point characterization [22] of the solution of the related optimization problem, which can lead to efficient algorithms. However, this discretization is difficult to extend to the CAUG due to its irregularly-distributed grid. We thus study the CAUG-based piecewise linear approximation of the anisotropic TV for reconstruction on the CAUG. We review the definition of the anisotropic TV [4,27]. Let D ⊂ R 2 be a bounded open set and W 1,1 (D) be the Sobolev space. For u ∈ W 1,1 (D), the anisotropic TV of u on D is defined by (17) u ATV := where | · | is the absolute value function, and (∇ x u, ∇ y u) is the weak gradient of u.
We further recall Courant's space [29] which is the space of continuous piecewise linear functions for the generated CAUG. For a continuous piecewise linear function v on Ω, we write the restriction of the function v on ∆ i ⊂ Ω as v| ∆i = v i and define (18) v We find that the corresponding gradient is given by ∇v i = [a i , b i ] , which is different from the gradient of a pixel depending on its adjacent pixels on the pixel grid. Based on the above discussions, we propose a regularizer suitable for the function v on Ω as the following form Motivated by the pixel-based piecewise constant approximation of the anisotropic TV, we rewrite R(v) as a composition of the 1 -norm and a matrix, leading to efficient algorithms. We next construct the matrix based on the grid G, and show that R(v) can be identified as the composition structure.
Proposition 5.1. Let v be the piecewise linear function defined by (18) for the grid G. If the regularizer R(v) is defined by (19), then there exists a matrix B ∈ R 2n×d such that the following equation holds where with respect to vertices on G.
Proof. We construct a matrix B by the definition of the anisotropic TV and Cramer's rule, and further prove that (20) holds as follows.
It can be verified that the anisotropic TV of each linear function is related to the gradient of the function and area of the corresponding grid element from definitions (17) and (18). For the function v i of the i-th grid element (i ∈ N n ), by definition (17) where int(∆ i ) is the largest open set that is contained in ∆ i . Let S i be area of the i-th grid element. Since a i and b i are the coefficients of the function v i with respect to (x, y) ∈ int(∆ i ), with area S i we have that (21) v i ATV = (|a i | + |b i |) · S i .
In particular, with three vertices of each grid element, the gradient of the linear function can uniquely determined by its values at these vertices through Cramer's rule. This leads to an underlying representation for the right-hand side of (21). We find that the gradient of each linear function can be represented by a matrixvector multiplication using Cramer's rule. Here the matrix is decided by the coordinates of vertices of each grid element and the vector is given by the values of the function at these vertices. By equation (8) we have that (x ik , y ik ) = T i (k), k ∈ N 3 , which are the coordinates of vertices of the i-th grid element. With the values of the function v i at these vertices, we obtain the gradient of v i through applying Cramer's rule to the linear equations In particular, the determinant of the coefficient matrix of the linear equations is equivalent to 2S i , and the gradient is given by (22) a Since equation (22) only involves the i-th grid element, we can rewrite it as a matrixvector multiplication in terms of the vector v. We then expand the 2 × 3 matrix to a 2 × d matrix by inserting zero vectors. Let C a i1 := y i2 − y i3 , C a i2 := y i3 − y i1 , (22), we define the 2 × d matrix as (23) B , Inverse Problems and Imaging Volume 14, No. 1 (2020),  where 0 c is the vector with all components equal to 0 for a positive integer c. We thus have the following formula for (22) as Moreover, by the 1 -norm the right-hand side of (21) has the following representation We construct the matrix B and prove that (20) holds using the above analysis. By (23) we can define the matrix B ∈ R 2n×d as (25) B which will be used to the following proof of (20). From (21) and (24), we have that and for each i ∈ N n it can be verified that Therefore, by (19) and the matrix B defined by (25) we find that In particular, for an image based on the pixel grid the above equation (20) reduces to the pixel-based piecewise constant approximation of the anisotropic TV if we choose B as the first-order difference matrix.
We now study the regularizer suitable for the solutionf described in Section 4. Sincef is represented by piecewise linear basis functions defined on G, we have the following corollary for the regularizer off by Proposition 5.1.
Corollary 5.2. Let B ∈ R 2n×d be the matrix defined by (25) for the grid G. If the solutionf of equation (9) is defined by (10) and the vector f ∈ R d is given by (14), then the regularizer R(f ) off can be identified as Proof. Since the solutionf can be represented by piecewise linear basis functions defined on G, we can write the restriction of the solutionf on ∆ i (i ∈ N n ) as f | ∆i =f i and definef i as (18). From (19), the regularizer off is denoted by Recalling B and f , we have that R(f ) = Bf 1 by Proposition 5.1.
We then formulate the regularized SPECT reconstruction problem as an optimization model via the penalized ML criterion [9,15] and the developed regularization method. The penalized ML estimate is achieved by maximizing the sum of the log-likelihood function of f , a negative regularization term via (26) and the indicator function of nonnegativity constraint [12]. Specially, we choose ϕ as the 1 -norm on R 2n , and compute the regularizer R(f ) (26) by the composition ϕ(Bf ) of the 1 -norm and the matrix B. We further define the logarithmic function of a vector x ∈ R q as ln x := [ln x i : i ∈ N q ] and the inner product of two vectors x, z ∈ R q as x, z := i∈Nq x i z i . Since the statistical model (16) guarantees a Poisson form of the underlying likelihood function, we develop an optimization model for solving (16) as the following formula In (27), 1 is the vector with all components equal to 1, β is a positive penalty weight, B ∈ R 2n×d is the matrix via (25), ϕ ∈ Γ 0 (R 2n ) (Γ 0 (R 2n ) denotes the space of all proper lower semi-continuous convex functions mapping from R 2n to R ∪ {+∞}, [1]), and ψ ∈ Γ 0 (R d ) given by is the indicator function of the nonnegativity constraint set R d + := {y ∈ R d : y ≥ 0}. The nonnegativity constraint is feasible and necessary since the radioactivity should be nonnegative at all vertices.
6. Iterative algorithm. We develop in this section an iterative algorithm for solving the resulting optimization model. We first characterize a solution of model (27) as a system of fixed-point equations. We then describe an iterative scheme to find the solution based on its fixed-point characterization. By the iterative scheme, we further develop an iterative algorithm which inherits the advantages of the fixedpoint proximity algorithm and the ML-EM algorithm.
For the description of the fixed-point characterization of a solution of model (27), we review some basic definitions and notation. Let S q + be the set of q × q symmetric positive definite matrices. For two vectors x, z ∈ R q and a matrix H ∈ S q + , the division is defined by x/z := [x i /z i : i ∈ N q , z i = 0], the weighted inner product is defined by x, z H := x, Hz and the corresponding norm is defined by x H := x, x 1/2 H . For a function ϑ ∈ Γ 0 (R q ), its proximity operator [12] with respect to a given matrix H ∈ S q + , denoted by prox ϑ,H , is a mapping from R q to itself, defined for a given vector x ∈ R q by prox ϑ,H (x) := argmin The subdifferential of ϑ ∈ Γ 0 (R q ) at a given vector x ∈ R q is the set defined by The subdifferential of the function ϑ and its proximity operator with respect to H ∈ S q + have the following relationship (28) Hz ∈ ∂ϑ(x) if and only if x = prox ϑ,H (x + z).
The conjugate ϑ * of ϑ ∈ Γ 0 (R q ) is defined at z ∈ R q by and an equivalent relationship of ∂ϑ and ∂ϑ * is given by (29) z ∈ ∂ϑ(x) if and only if x ∈ ∂ϑ * (z) for x ∈ dom(ϑ) and z ∈ dom(ϑ * ). For a discussion of these relations, see, e.g., [1] or [14]. Therefore, a solution of model (27) may be characterized as a system of fixed-point equations in the following theorem using a way similar to [12,14]. (27), then for any Q ∈ S d + and P ∈ S 2n + , there exists a vector h ∈ R 2n such that the following coupled equations hold Conversely, if there exist Q ∈ S d + , P ∈ S 2n + , f ∈ R d and h ∈ R 2n satisfying the coupled equations (30)- (31), then f is a solution of model (27).
Conversely, if there exist Q ∈ S d + and P ∈ S 2n + such that (f , h) ∈ R d+2n satisfies (30)- (31), then all the arguments discussed above are reversible.
We rewrite two equations in Theorem 6.1 as a compact formula for development and convergence analysis of the underlying iterative scheme. To this end, an operator T : R d × R 2n → R d × R 2n at a vector w := (f , h) ∈ R d × R 2n , used to integrate the two proximity operators in (30)- (31), is defined by Furthermore, we define an operator R : where ∇φ(f ) is the gradient of the function the compact formula of the coupled equations (30)- (31) can be denoted by where I d and I 2n are two identity matrices. In particular, let Ψ(w) := λψ(f ) + (λβϕ) * (h) and (38) D := diag(Q, P ), this implies that the operator T is considered as the proximity operator of the function Ψ with respect to the matrix D, denoted by T = prox Ψ,D . The operator T is firmly nonexpansive with respect to D by [14], which means that . Moreover, we have the following property in terms of matrices G and D.
Lemma 6.2. If the matrix G is defined by (36) and the matrix D is defined by (38), then G D > 1.
Proof. Recalling matrices G, D and following the proof of Lemma 3.2 in [14], for w := (f , h) ∈ R d × R 2n we have that This leads to the desired result.
Lemma 6.2 shows that the matrix G is expansive. As indicated in [22], a simple Picard iteration based on the fixed-point equation (37) may not be convergent due to the expansiveness of the matrix G. We then describe an implicit iterative scheme based on the fixed-point characterization (37). In particular, the underlying implicit iterative scheme can be implemented explicitly and guaranteed to converge under proper conditions. To this end, let we have the following splitting strategy for the matrix G, denoted by and observe that the fixed-point equation (37) can be rewritten as This leads to an implicit iterative scheme for finding a fixed point as the following formula which is described in detail by where N 0 := {0, 1, 2, . . .}. The resulting iterative scheme overcomes the difficulty arising from the non-smoothness of model (27). Moreover, it involves the preconditioning strategy, aiming at obtaining a numerical solution in less time via the preconditioning matrix applied for matrix computing.
We now develop an effective iterative algorithm from the resulting iterative scheme through selecting a proper preconditioning matrix. To this end, we introduce a nontrivial selection of the preconditioning matrix Q motivated by the EM algorithm [12,13]. Indeed, we choose P := I 2n and A ij for the matrix A, γ > 0, and ω k j is a positive parameter varying with iteration number k ∈ N 0 . We may set ω k j as the current vertex estimate f j at the first few iterations, and fix it after certain iteration. This gives the following preconditioned fixed-point proximity algorithm.

Algorithm 1 Preconditioned Fixed-point Proximity Algorithm (PFP 2 A)
The above algorithm can be applied to reconstruction on the CAUG or on the pixel grid by choosing different matrices A and B. For implementation of the resulting algorithm, by [12,22] we have that prox λψ,Q (f ) j = max{f j , 0}, j ∈ N d for f ∈ R d , and 7. Convergence analysis. We provide in this section convergence analysis for the PFP 2 A motivated by [17,32]. For convenience in the convergence analysis of the resulting algorithm, we first transform the implicit iterative scheme (41) to an explicit iteration through defining an operator [14]. Suppose that for any u ∈ R d+2n , there exists a unique w ∈ R d+2n such that We may define the operator F : R d+2n → R d+2n as where F has no closed form. This means that the implicit iterative scheme can be symbolically represented as the following explicit iteration With the explicit iteration (44), we will apply the following convergence result of an explicit iterative scheme to prove convergence of the resulting algorithm.
Proposition 7.1. Let U be the fixed-point set of the operator L : R d+2n → R d+2n , {v k ∈ R d+2n : k ∈ N 0 } be the sequence generated by the iterative scheme v k+1 = L(v k ) for any initial value v 0 ∈ R d+2n , and W 1 , W 2 ∈ S d+2n + . If the following conditions hold (i) the operator L is continuous, , then the sequence {v k : k ∈ N 0 } converges to a fixed point of L.
Proof. We first prove that the sequence {v k : k ∈ N 0 } converges a cluster point via condition (ii). Using condition (ii), we find that for any There exists a subsequence {v kp : p ∈ N 0 } that converges to a cluster point v * ∈ R d+2n . Moreover, we sum the inequality in (ii) for k from 0 to N − 1, and have that which means that lim k→+∞ v k − v k+1 W2 = 0. Therefore, the sequence {v k : k ∈ N 0 } has at most one cluster point.
We then prove that the cluster point is a fixed point of the operator L by condition (i). Since the operator L is continuous and v k+1 = L(v k ), the cluster point v * is a fixed point of L and v * ∈ U.
As a result, the sequence {v k : k ∈ N 0 } converges to a fixed point of L.
We need to introduce two lemmas before applying Proposition 7.1 to prove the convergence of the PFP 2 A. The first lemma describes the Lipschitz continuity of the gradient of φ using a similar way to [12]. For two vectors x, z ∈ R q , the componentwise multiplication is defined by Lemma 7.2. For A ∈ R m×d , g ∈ R m and r ∈ R m in (16), if φ is the function on R d defined by (35), then ∇φ is Lipschitz continuous with the constant η, where

It can be verified that
This proves the desired result. We then introduce the second lemma for analyzing the symmetric positive definiteness of the related matrices arising from the PFP 2 A. Notation is required to analyze that. For the selected matrix Q, we set ω k j := ω j (j ∈ N d ) for k > k 0 which means that we fix the variable parameter after k 0 iterations for the PFP 2 A, where ω j > 0. For the selected matrices Q and P = I 2n in the PFP 2 A, by the defined matrices D (38) and G 1 (40) we have that Let W := DG 1 and (45) W := W − L with L := diag(2ληI d , 0) for two parameters λ, η, which will be used later. Let Based on the above notation, we have the following lemma about two matrices W and W . Proof. By definition (45), the symmetric positive definiteness of W implies that of W . It is obvious that W is a symmetric matrix. We need to prove that W is positive definite. We first prove that the matrix Q − 2ληI d is positive definite for the following proof of W . Since we find that all eigenvalues of the matrix I d − 2ληQ −1 are not less than ε by the second condition in (47). This means that Q − 2ληI d is a positive definite matrix. We then prove that W is positive definite. Let It can be verified that for W , that is, W and diag(I d , I 2n − EE ) are congruent. This implies that W is positive definite if and only if E 2 < 1. We now verify that E 2 < 1. Since 2λη/γ ≤ (1 − ε)ξ, it can be verified that for the matrix Q − 2ληI d , its each diagonal element satisfies This means that (Q − 2ληI d ) −1/2 2 ≤ (εξγ) −1/2 . Hence, by the first condition in (47) we observe that This proves the desired results.
Convergence analysis is performed for the PFP 2 A based on the above analyses. By Proposition 7.1, the resulting algorithm is convergent if we can prove that the operator F is continuous and an underlying inequality holds. Here we employ the closed graph theorem to prove the continuity of the operator F, and let gra(F) denote its graph defined by gra(F) := {(u, w) : u, w ∈ R d+2n , w = Fu}. The convergence analysis of the PFP 2 A can be given by the following theorem.
Theorem 7.4. Let {w k = (f k , h k ) : k ∈ N 0 } be the sequence yielded by the PFP 2 A, and U be the fixed-point set of the operator F defined by (43). If parameters λ and γ satisfy (47), then the sequence {w k : k ∈ N 0 } converges to a fixed point of the operator F and {f k : k ∈ N 0 } converges to a solution of model (27).
Moreover, it can be verified that . Recalling definition (45) of W , by Lemma 7.3 it is a symmetric positive definite matrix. We then obtain that Therefore, the sequence {w k : k ∈ N 0 } converges to a fixed point of F via Proposition 7.1. Moreover, the fixed point is also a solution of the coupled equations (30)-(31) by the definition of F, and hence, {f k : k ∈ N 0 } converges to a solution of model (27). 8. Numerical experiments. We perform in this section numerical experiments to demonstrate the effectiveness of the developed method for 2D SPECT reconstruction. We compare the developed method with the ML-EM [7], the MAP-EM [3] on the CAUG, and the discrete TV regularized reconstruction method on the pixel grid in terms of the reconstructed images quality and computational performance. In particular, for image visualization the related radioactivity distributions are colored via ImageJ tool.
We apply the chest phantom with size 128 × 128 to yield projection data with different noise levels. We obtain noise-free parallel-beam SPECT data (figure 2(c)) with 120 angles and 185 detector bins from the phantom (figures 2(a)-(b)) through the matrix-vector multiplication on the pixel grid via MATLAB. Based on the simulated data, we use a Poisson random number generator to further create projection data with low (8%) and high (11%) Poisson noise ( figure 3(a)). With the given attenuation distribution ( figure 2(b)), we obtain the system matrix A through calculating the formula (12) via Simpson's rule, which is similar to [21]. We set the vector r with all components equal to 10 for background activity.  We then generate CAUGs used to SPECT reconstruction. We seek the initial estimates with important features. These estimates are reconstructed by the discrete TV regularized reconstruction method on the pixel grid with 200 iterations through C++. In particular, the obtained initial estimates are regarded as the initial values used to reconstruction on the CAUG. Here the discrete TV regularized reconstruction method on the pixel grid is the PFP 2 A with the pixel-based piecewise constant approximation of the anisotropic TV (DTV). In the context of the pixel grid, A is  matrices I 128 and B. We further generate the corresponding CAUGs (figure 3(c)) by the method described in Section 3 through MATLAB. It takes about 3.6 seconds to generate a CAUG. Here we choose different 1 from 10 −4 to 10 −3 according to different initial estimates, set 2 = 0.1 and k = 0 in (7). Due to the use of the given attenuation distribution, there are no obvious differences between CAUGs for three initial estimates. The use of the CAUG substantially reduces the dimension of the solution space (table 1)    The results (figures 4-5) are reconstructed from projection data with different noise levels through four competing methods using C++. These methods include the developed method, the ML-EM and the MAP-EM on the CAUG, and the discrete TV regularized reconstruction method on the pixel grid. The developed method is the PFP 2 A with the proposed regularization method defined on the CAUG (RUG). When using the PFP 2 A, we fix the preconditioning matrix after 200 iterations, and set ω k j (j ∈ N d ) as the estimate f 200 j when k > 200. In particular, the step-size parameter λ/γ in the PFP 2 A is very small due to the huge Lipschitz constant, this results in very slow convergence. By the EM-preconditioner we may choose the parameters γ > B 2 2 /(εξ) and λ > 0 such that λ/γ is from 1 to 2. We further choose the regularization parameter from 20 to 800 such that the competing methods yield the best normalized mean square error (NMSE) values. Iterations for reconstruction stop when the relative error of the iterative sequence {f k : k ∈ N 0 } is less than 10 −5 , given by f k+1 − f k 2 / f k+1 2 < 10 −5 . For image visualization, the discrete forms of the reconstructed functions on the CAUG are their samplings on the pixel grid with size 128 × 128.
We now evaluate quantitatively the reconstructed images using some metrics due to nonlinear and nonstationary noise property of the SPECT system. To this end, we select the reconstructed images from projection data with high noise as the assessed objects. We further select the region of interest (ROI) and the line profile in the reference image ( figure 7(a)), and the following metrics (figure (a) (b) Figure 6. (a) The relative error of the iterative sequence varies with computing time (starting with 0.6 seconds) for reconstruction on the CAUG and on the pixel grid from projection data with high noise; (b) Total time of obtaining the reconstructed images, including 6.8 seconds for yielding the initial estimate and 3.6 seconds for generating the CAUG.
7(b)-(f)) suitable for such a system. The NMSE is defined by NMSE := 100 × with Z as the reconstructed image and X as the reference image. For the given ROI, the signal-to-noise ratio (SNR) is defined by SNR := ξ ROI /σ ROI , where ξ ROI is the averaged activity and σ ROI is the standard deviation in the ROI. The contrast-to-noise ratio (CNR) between two ROIs is defined by CNR := |ξ ROI1 − ξ ROI2 |/ σ 2 ROI1 + σ 2 ROI2 . The full width at half maximum (FWHM) is defined by FWHM := 2 √ 2 ln 2σ, where σ is the standard deviation of the fitted Gaussian function for the line profile. Finally, we consider computational costs of the four competing methods (figure 6).
We further conduct numerical experiments for an underdetermined system based on the pixel grid, compared to the above overdetermined systems for SPECT imaging. We first yield projections with 120 angles and 128 detector bins using the same way as in the above experiments. In particular, the resulting projections may be considered as the truncation cases of the above projections with size 120 × 185 through removing its first 29 and last 28 columns. We then obtain the reconstructed images (figure 8) from the resulting projections with different noise levels through the DTV and the RUG. Moreover, NMSE, CNR and FWHM (table 2) are used to evaluate the reconstructed images from the resulting projection with high noise.  The above results demonstrate the superiority of our method over the other competing methods in improving image quality and reducing computational costs. In particular, by table 1 an image can be represented by using far fewer CAUGbased basis functions than pixel-based basis functions. Visual comparisons from figures 4-5 and 8 show that the developed method can effectively suppress noise and preserve edges of the reconstructed images. From figure 6, we see that the regularized reconstruction method on the CAUG efficiently reduces computational costs. Quantitatively evaluations in figure 7 and table 2 further show that our method outperforms the other competing methods in terms of these categories for the reconstructed results from projection data with high noise. 9. Conclusions. We developed a novel SPECT reconstruction method in this paper. It consists of the CAUG-based piecewise linear approximation of an integral equation model for SPECT, the edge-preserving regularization method defined on the CAUG, and the preconditioned fixed-point proximity algorithm to solve the resulting optimization problem. Numerical experiments demonstrate that the developed method outperforms the other competing methods in terms of NMSE, noise suppression (standard deviation, SNR and CNR), and local spatial resolution (FWHM). In particular, the CAUG-based reconstruction method can significantly reduce computational costs, as compared to the reconstruction method on the pixel grid.