Variational regularization of the weighted conical Radon transform

Recovering a function from integrals over conical surfaces recently got significant interest. It is relevant for emission tomography with Compton cameras and other imaging applications. In this paper, we consider the weighted conical Radon transform with vertices on the sphere. Opposed to previous works on conical Radon transform, we allow a general weight depending on the distance of the integration point from the vertex. As first main result, we show uniqueness of inversion for that transform. To stably invert the weighted conical Radon transform, we use general convex variational regularization. We present numerical minimization schemes based on the Chambolle-Pock primal dual algorithm. Within this framework, we compare various regularization terms, including non-negativity constraints, $H^1$-regularization and total variation regularization. Compared to standard quadratic Tikhonov regularization, TV-regularization is demonstrated to significantly increase the reconstruction quality from conical Radon data.


Introduction
In this paper, we study the problem of reconstructing a function from its weighted conical Radon transform gf@z; A a Z S@z; A f@xA U @kx zkA dS@xA : Here S@z; A R n is a conical surface with vertex z on the unit sphere S n I , semi-axis z and half-opening angle P H; =P, U is a radial weight function, and dS the standard surface measure.
There are various variants of conical Radon transforms that differ in the restrictions of the opening angle, the axis direction or the locations of the vertices of the conical surfaces, and might include different weights. Many of such versions of conical Radon transforms have been studied in the literature, see [1][2][3][4][5][6][7][8][9][10][11][12][13][14][15][16][17][18]. For the special case of a weight U@rA a r m with m P Z, our transform gf has been studied in [19] in general dimension, and in [20] for the weight U@rA a e r in dimension n a P. Some other works (for example [3,12,17]) consider weights of the form U@rA a r m for overdetermined conical Radon transforms. The ill-posed character is microlocally analyzed in [21] in dimension n a P for the case U@rA a I, where also artifacts are characterized. The use of a general weight for the conical Radon transform has not been previously studied for any variant of the conical Radon transform.
The inversion of conical Radon transforms arises in single photon emission tomography (SPECT) using Compton cameras [22,23], as well as in other applications such as single scattering optical tomography [24]. In the context of SPECT with Compton cameras, the general radial weight that we consider in this paper allows to include the physically relevant effect of attenuation of photons. In the context of SPECT using the classical Anger cameras where the exponential and the attenuated Radon transform appear this issue is well investigated; see for example [25][26][27][28][29][30][31][32][33][34]. For the conical Radon transform, attenuation is investigated much less and we are only aware of the works [35,36] considering the special case U@rA a e r .
The main contributions of this article are as follows. First, we derive a uniqueness result for g extending the result of [19] to general weights. Second, in order to address the ill-posedness of inversion of the weighted conical Radon transform, we apply convex variational regularization g; @fA Xa I P kgf gk P L 2 C ©@fA 3 min f : Here denotes a positive regularization parameter and © is a convex regularizer that allows to include positivity as well as total variation (TV) regularization. For the implementation of the numerical minimization, we apply the primal dual algorithm of Chambolle and Pock [37]. To the best of our knowledge, neither general convex variational regularization nor the Chambolle-Pock algorithm have previously been studied for the conical Radon transform. We present numerical studies comparing various regularizers. The presented results clearly suggest that TV-regularization outperforms the other regularizers for the considered phantom type.

Outline
This paper is organized as follows. In Section 2, we define the weighted conical Radon transform and study its mathematical properties. In Section 3, we consider its variational regularization. Additionally, we derive the numerical minimization algorithm and present numerical results for L P -regularization, H I -regularization and TV-regularization, all with and without positivity constraint. The paper concludes with a short summary and outlook presented in Section 4.

The conical Radon transform
Throughout this paper, U P C I @H; IA; RA denotes a given function with U@rA ! H for all r ! H. For z P S n I and P H; =P, we denote by S@z; A a fz C r! j r ! H and ! P S n I with ! z a os g (2.1) the surface of a right circular half cone in R n with vertex z, central axis z, and half opening angle . We write D X Xa C I H @B I @HAA for the space of smooth functions with compact support in the ball B I @HA, and D Y Xa C I @S n I ¢ H; =PA. The results in this section extend the results of [19] for the conical Radon transform with weights U@rA a r m to the case of general radial weights.

Definition and basic properties
We first define the transform studied in this paper. The weighted conical Radon transform integrates a function f P D X supported inside the unit ball over cones with vertices on the unit sphere S n I a fx P R n j kxk a Ig and central axis orthogonal to S n I pointing to the interior of the sphere. Alternative expressions for the conical Radon transform that will be helpful for our analysis are derived next.

Adjoint transform
For g P D Y we define the weighted conical backprojection g ] g X R n 3 R by g ] g@xA Xa Z S n 1 Z =P H g@z; A U @kx zkA ¢ @@x zA z C kx zk os@ AA d dS@zA ; (2.5) if x P B R0 @HA, and g ] g@xA a H otherwise. In the next proposition, we will show that g ] is the (formal) adjoint of g with respect to the L P -inner products hf I ; f P i L 2 Xa Z R n f I @xAf P @xA dx where we used the Theorem of Fubini.

Decomposition in one-dimensional integral equations
Next we derive an explicit decomposition of the conical Radon transform in one-dimensional integral operators (see Theorem 2.4). For that purpose we use the spherical harmonic decompositions f@rA a I X aH N@n;`A X kaI f`; k @rAY`; k @A ; (2.8) @gfA@z; A a I X aH N@n;`A X kaI @gfA`; k @ AY`; k @zA : Here Y`; k , for`P N and k P fI; : : : ; N@n;`Ag, denote spherical harmonics [38,39] of degree`forming a complete orthonormal system in S n I . The set of all @`; kA with`P N and k P fI; : : : ; N@n;`Ag will be denoted by I@nA. Let C denote the Gegenbauer polynomials normalized in such a way that C @IA a I. As in [19], we derive different relations between f`; k and @gfA`; k in the form of Abelian integral equations.
Theorem 2.4 (Generalized Abel equation for f`; k ). Let f P D X and let f`; k and @gfA`; k be as (2.8) and (2.9) for @`; kA P I@nA. Then, for P H; =P, @gfA`; k @ A a jS n P j Z I sin@ A f`; k @A K`@ ; A q P @sin@ AA P d ; @sin@ AA n I @sin@AA n P @sin@ C AA n C @n PA=P @os@AAd : (2.12) Splitting the integral in one integral over < =P and one over ! =P and proceeding similar as in [19, Theorem 3.2] yields the claim.
Proposition 2.5. Let f P D X and let f`; k and @gfA`; k be as (2.8) and (2.9).
Further, for every @`; kA P I@nA denote (a) g`; k @tA Xa jS n P j I @I tA @n PA=P @gfA`; k @ros p tA;

Uniqueness of reconstruction
In this section, we show uniqueness of recovering the function by its conical radon transform and thus the injectivity of g by showing solution uniqueness of the Abelian integral equations in Theorem 2.7. Since the kernel function of the Abelian integral equation (2.10) has zeros on the diagonal, the proof of the uniqueness relies on the uniqueness result derived in [19], which we briefly state at this point: For a; b P R with a < b we set ¡@a; bA Xa f@t; sA P R P j a s t bg.
Lemma 2.6. Suppose that F X ¡@a; bA 3 R, where a < b, satisfies the following: (F1) F P C Q @¡@a; bAA. (F2) N F Xa fs P a; bA j F @s; sA a Hg is finite and consists of simple roots.
(F3) For every s P N F , the gradient @ I ; P A Xa rF@s; sA satisfies I C I P I I C P > H: (2.14) Then, for any g P C@a; bA, the integral equation Vt P a; bX R t a F @t;sA p t s f@sAds a g@tA has at most one solution f P C@a; bA. For any f P D X and any @`; kA P I@nA, the spherical harmonic coefficient f`; k of f can be recovered as the unique solution of V P H; =PX @gfA`; k @ A a jS n P j Z I sin@ A f`; k @A K`@ ; Ad q P @sin@ AA P ; with the kernel functions K`defined by (2.11).
Proof. We proceed as in the proof of Theorem 3.5 in [19]. Let f P D X vanish outside a ball of radius I a P , where a P @H; IA. We show that equation (2.13) has a unique solution, which is sufficient according to Lemma 2.5. For that purpose, we verify that F`X ¡@a; IA 3 R satisfies the conditions (F1)-(F3) in Lemma 2.6.
Ad (F1): Using the abbreviations U n for the function r U 3 U@rA r n P and C Xa C @n PA=P , the kernel F`can be written in the form The function F`is clearly smooth on f@t; sA P ¡@a; IA j t T a sg. Using C@ xA a @ IA`C@xA and that U n is analytic, F`is an even real analytic function in p t s which shows that F`is also smooth on the diagonal f@t; sA P ¡@a; IA j t a sg.
Ad (F2): Let v@sA Xa F`@s; sA a PU n @ p sAC@ p I sA denote the restriction of the kernel to the diagonal.
As an orthogonal polynomial, the function C has only a finite number of isolated and simple roots. We conclude the same holds true for the function v, using that U is a positive function.

Convex regularization and iterative minimization
The stability analysis made in [20] (performed there in 2D) shows that the inversion of the conical Radon transform is highly ill-conditioned (compare [21], where the artifacts are characterized in 2D). To address the ill-posedness, in our previous works [19,20] we used standard quadratic Tikhonov regularization. In this paper, we go one step further and apply variational regularization allowing a general convex regularization term. For the following we set X Xa L P @B I @HAA and Y Xa L P @S n I ¢ H; =PA.

Variational regularization
In order to stably solve gf a g, we consider variational regularization which consists in minimizing the generalized Tikhonov functional g; @fA Xa I P kgf gk P L 2 C ©@fA : Here is some non-negative constant and ©X X 3 H;I is a proper convex, coercive and weakly lower semi-continuous functional. Tikhonov regularization is used to stably approximate ©-minimizing solutions of g@fA a g, which are defined as elements in rg minf©@fA j f P X g@fA a gg.
We first derive the boundedness of f U 3 gf with respect to the L P -norm.
Lemma 3.1. Suppose r U 3 r n Q U@rA P is integrable over @H;PA. For some constant c P @H; IA we have kgfk L 2 ckfk L 2 for all f P D X . In particular, g can be uniquely extended to a bounded operator g X X 3 Y, which has bounded adjoint.
Proof. Let z P S n I and Q P O@nA satisfy Qe I a z. Integration over z P S n I and using that R P H r n Q U@rA P dr is finite yields the claimed estimate. Because D X X is dense, the operator can be extended to a bounded operator on X in a unique way. The adjoint is therefore bounded, too.
The continuity of g together with standard results for convex variational regularization implies the following.
Theorem 3.2 (Existence, stability and convergence). Let ©X X 3 H; I be proper, convex, coercive and weakly lower semi-continuous, and suppose r U 3 r n Q U@rA P is integrable over @H;PA. Then the following hold: (a) For every g P Y and > H,¨g ; has at least one minimizer. (b) Let > H, g P Y, and let @g k A kPN P Y N satisfy kg g k k L 2 3 H.
Then, every sequence f k P rg min¨g k ; has a weakly convergent subsequence @f @kA A kPN , and its limit f is a minimizer of¨g ; with ©@f @kA A 3 ©@fA for k 3 I. (c) Let g P rn@gA, @ k A kPN P @H;IA N converge to zero, and let @g k A kPN P Y N satisfy kg g k k k . Suppose further that @ k A kPN P @H;IA N satisfies k 3 H and P k = k 3 H as k 3 I, and choose f k P rg min¨g k ; k for every k P N. @f k A kPN has a weakly converging subsequence. The limit f of every weakly convergent subsequence @f @kA A`P N of @f k A kPN is a ©-minimizing solution of g@fA a g and satisfies ©@f @kA A 3 ©@fA. If the ©-minimizing solution of g@fA a g is unique, then f k * f.
Proof. Follows from the boundedness g derived in Proposition 3.1 together with general results for convex variational regularization [40].
Note that the uniqueness of a ©-minimizing solution of g@fA a g is guaranteed if its solution is unique. Uniqueness of the ©-minimizing solution is also guaranteed in the case where © is strictly convex. This property is satisfied, for example, for standard Tikhonov regularization where ©@fA a kfk P L 2 , or for q -regularization where ©@fA a P P£ jhf ; ' ij q with q > I and some frame @' A of X.

Iterative minimization
In the numerical results presented below, we consider the following instances for the regularizer ©X X 3 H; I: ©@fA a I P kfk P L 2 (L P -regularization); ©@fA a I P kDfk P L 2 (H I -regularization); ©@fA a kDfk L 1 (TV-regularization).
Additionally, we consider any of these methods with an added convex constraint. All resulting regularization functionals are proper, convex, coercive and weakly lower semi-continuous.
For numerically minimizing the Tikhonov functional (3.1), we consider its discrete counterparts. For that purpose, let f P @R NCI A n be the discrete phantom, g P R P ¢@QCIA discrete data and C X @R NCI A n 3 R P ¢@QCIA the discretization of the forward operator. For all considered regularizers, we can write the resulting discrete Tikhonov functional in the form Φ@f A Xa I P kCf g k P P C q kLfk q q C I M @fA : Here I M denotes the indicator of some convex set M @R NCI A n defined by I M @fA a H if f P M and I M @fA a I else. In particular, in the case that we take M a @H;IA NCI A n it guarantees non-negativity. The mapping L P fD; Ig stands either for the discrete gradient D X @R NCI A n 3 @@R NCI A n A P in the case of H I -regularization and TV-regularization, or for the identity operator I on @R NCI A n in the case of L P -regularization. The parameter q P fI; Pg is taken q a I in the case of TV regularization, and q a P in the cases of L P -regularization or H I -regularization. p kCI 2 @p k C @Cu k g AA=@I C A 6: q kCI 2 @q k C Du k A=mxf1;jq k C Du k jg 7: f kCI 2 P M @f k C £ p kCI D £ q kCI A 8: u kCI 2 f kCI C @f kCI f k A 9: k 2 k C I

10: end while
The Tikhonov functional (3.2) can be minimized by various convex optimization methods [41]. In this work we use the minimization algorithm of [42], which is a special instance of the Chambolle-Pock primal-dual algorithm [37]. It can be applied to any instance of regularization functional that we consider in this paper. For TV-minimization with convex constraint, the resulting al- Left: The numerical phantom f P f P R PSU¢PSU consists of the superposition of ellipse and two star shaped tumor like structures. Right: Simulated data g a Cf C ξ P R PHH¢ISI with S7 noise added.
and H I -regularization, the algorithms look similar, only line 6 has to be replaced by the update rule q kCI 2 @q k C Lu k A C : For motivation and a derivation of Algorithm 1 as well as a convergence analysis we refer to the original papers [37,42].

Numerical simulations
For the presented numerical results, we assume weight U@rA a e r and consider the two dimensional case. In this case, the weighted conical Radon transform reduces to the attenuated V-line transform @gfA@'; A a X a¦I Z I H f @@os@'A; sin@'AA r@os@' A; sin@' AAA e r dr ; (3.3) where ' P H;PA and P H; =P. The V -line transform (3.3) corresponds to transforms studied in [19] and [17].

Discrete forward and adjoint operator
The discrete forward and adjoint operators are defined by considering the entries f i I ; i P as sampled values of a function f at locations @ I; IA C P@i I ; i P A=N for @i I ; i P A P fH; : : : ; Ng P and replacing the function f by a bilinear interpolant T f. The attenuated V-line transform Cf a g is discretized by numerically computing the integral of the interpolant T f over each of the two branches of the V-lines with the composite trapezoidal rule. We take @os@'kA;sin@'kAA with 'k a P@k IA=P for k P fI; : : : ; P g for the vertex positions, and ` a `=@PQA with`P fH; : : : ; Qg for the opening angles. For the numerical integration of each branch, we use N CI equidistant radii in the interval H;P. The discrete backprojection operator is defined by using the trapezoidal rule similar to the well-known approach for the classical Radon transform [43].
For the following numerical results we use N a PST, P a PHH, Q a ISH, and a H:S. The used discrete phantom f and the numerically computed data Cf with added Gaussian noise are shown in Figure 3.1.

Reconstruction from exact data
We first investigate the variational regularization methods on simulated data without noise added. We consider L P -regularization, H I -regularization and TV-regularization. All methods are used with and without positivity constraint. For comparison purposes, we also perform reconstructions using plain least squares, positivity constrained least squares and the direct Fourier reconstruction method from [20]. The regularization parameter has been taken a H:HHP for TV and H I regularization and a H:HI for L P regularization.
The reconstruction results after 700 iterations are shown in Figure 3.2. Total variation minimization (with and without positivity constraint) clearly outperforms all other methods. In particular, L P -as well as H I -regularization contain a ghost source corresponding to the big disc, which is not contained in the TV-reconstruction. Using the least squares methods, the ghost image is also not that severe.
To quantitatively evaluate the reconstructions, we compute the relative squared P -error and relative squared`P-residuals, respectively, E P @f j A Xa P i1;i2 jf j i I ; i P f i I ; i P j P P i1;i2 jfi I ; i P j P (3.4) R P @f j A Xa P k;`j Cf j k;` gk;`j P P k;`j gk;`j P : (3.5) Logarithmic plots of E P @f j A and R P @f j A are shown in Figure 3.3. The reconstruction error for the TV-regularization is much smaller than for all other methods. The residuals, as expected, are smallest for plain least squares. As a consequence of the ill-posedness, this does not imply small reconstruction error. In fact, the least squares methods show a slight semi-convergence behavior which is due data error introduced by the numerical implementation and the ill-posedness of the inverting the attenuated V-line transform.

Reconstruction from noisy data
To test the algorithms in a more realistic situation, we repeated the simulation studies where we added additive Gaussian noise ξ P R @P CIA¢Q to the simulated    Logarithm log IH E P @f j A of squared reconstruction error. Right: Logarithm log IH R P @f j A of squared residuals.
data. The relative`P-error a q P k;`j ξk;`j P q P k;`j gk;`j P is approximately S7. We again consider L P -regularization, H I -regularization, TV-regularization, least-squares minimization and the Fourier reconstruction method from [20], where each method is used with and without positivity constraint. In order to stabilize the reconstruction, the regularization parameter had to be increased for all methods. It has been taken a H:IR for L P -regularization, a H:HT for H I -regularization, and a H:HIS for TVregularization.
Reconstruction results from noisy data are shown in Figure 3.4. We used 200 iterations for the regularized iterations and stopped the least squares iterations after 15 iterations in order to prevent overfitting. Logarithmic plots of the squared relative`P-error E P @f j A and the squared relative`P-residual R P @f j A are shown in Figure 3.5. As can be seen from the reconstruction as well as the evolution of the reconstruction error, TV-regularization again performs best for the considered type of phantom. In particular, in all other methods, the ghost image is clearly visible and the reconstruction error is much larger than for TV-regularization.

Conclusion and outlook
In this paper, we investigated variational regularization methods for the stable inversion of the conical Radon transform with vertices on the sphere. We presented convergence results of variational regularization and a numerical minimization algorithm based on the Chambolle-Pock primal dual algorithm [37,42]. In particular, the algorithm framework allows a TV-regularizer,    Logarithm log IH E P @f j A of squared reconstruction error. Right: Logarithm log IH R P @f j A of squared residuals.
quadratic penalties, positivity constraint as well as their combinations. In terms of relative`P-reconstruction error, we found TV-regularization to clearly outperform least squares quadratic-regularization (all with or without positivity constraint).
The results in this paper show that compared to quadratic regularization, the use of non-smooth convex regularization can clearly improve reconstruction of the conical Radon transforms. Several interesting generalizations of variational regularization are possible. First, one could replace the sphere with more general surfaces of vertices and allow different axis for each vertex location. Second, the algorithmic framework can be extended to variable axis on the sphere. The transform then depends on the Pn I independent variables and is therefore highly overdetermined. In applications such as emission tomography with Compton camera, the noise level is usually high and the full overdetermined transform must actually be used for image reconstruction. Moreover, noise statistics as well as factors governing system performance such as Doppler broadening or polarization [22] should be included in the forward and adjoint problem. Resulting algorithms based on conical Radon transforms (which require binning of Compton camera raw data) should be compared with list mode EM reconstruction algorithms [44] that work without data binning.
Besides numerical investigations, several theoretical aspects have to be addressed in future work. This includes stability and uniqueness analysis as well as range descriptions of the various forms of conical Radon transforms.