Filtered backprojection (FBP) reconstruction algorithms for parallel-beam projection data with and without attenuation have been well established in [1−6]. Without attenuation, the FBP algorithm has been extended to fan-beam data acquisition geometry in [7, 8], and has been widely used in Computed Tomography (CT). In Single Photon Emission Computed Tomography (SPECT), the attenuation needs to be compensated for in quantitative studies. However, the existence of an FBP type algorithm for attenuated fan-beam projections had been an open problem for a while. In [9], the authors presented a circular harmonic decomposition (CHD) reconstruction algorithm for SPECT with fan-beam collimation and uniform attenuation compensation. The CHD method does not indicate any knowledge about the existence of an FBP algorithm for fan-beam projections. Recently, an explicit inversion formula of FBP type was obtained for attenuated fan-beam projections in [10], but the derivation seems to be quite complicated to people in the area of engineering and computer science. In this paper, based on Novikov’s inversion formula and a partial derivative relation between Cartesian and polar coordinates, a FBP reconstruction algorithm for attenuated fan-beam projections is obtained. The proof seems to be quite elementary compared with other existing work. The algorithm can be reduced to the special case of uniform attenuation since Novikov’s inversion formula is valid when the attenuation map is constant in a unit disk [5] or is zero.

This paper is organized as follows. Section II reviews Novikov’s inversion formula for parallel-beam geometry. Section III describes the derivation of a fan-beam FBP reconstruction algorithm that compensates for non-uniform attenuation. Discretization of the FBP algorithm and computer simulations are presented in Section IV. Section V concludes the paper with discussions and comments.

Let R² denote the two-dimensional (2D) Euclidean space with point representation (x, y) or in the Cartesian coordinates and (r, ) in the polar coordinates, respectively. For each function f (x, y) in R², f_θ(s, t) = f_θ(s cosθ − t sinθ, s sinθ + t cosθ) denotes the same function in the Cartesian coordinates (s, t)_θ that is the rotation of (x, y) by an angle θ along the counter-clockwise direction. The relation between these two coordinates can be depicted by the following transformations

(1)

The coordinate (x, y) is defined in the object space, while the coordinate (s, t)_θ is associated with the rotating detector. Hereafter, we call (x, y) as object coordinate, (s, t)_θ as parallel-beam detector coordinate, and (s, θ) or (l, θ) as the parallel-beam projection coordinate. For simplicity, we introduce two vector symbols of

(2)

and one complex symbol of .

In SPECT imaging, f (x, y) represents the distribution of radiotracer concentration inside the body. The γ-ray photons, emitted from the radiotracer, are attenuated inside the body by a linear attenuation coefficient distribution μ(x, y) before they arrive at the detector. Let p(s, θ) be the accumulated photon counts at position s and detector view angle θ in the parallel-beam geometry, then

(3)

where . Function p(s, θ) is the so-called attenuated Radon transform, and will be denoted by R_μf. It has been known that f (x, y) can be exactly reconstructed from p(s, θ) using Novikov’s inversion formula in parallel-beam geometry. Before presenting that formula, we introduce some mathematical definitions below. Define the Hilbert transform as

(4)

the Radon transform as

(5)

and the derivative of the Hilbert transform as

(6)

Using the conventions of [4], Novikov’s inversion formula can be expressed as

(7)

where , and l indicates the rays in the parallel-beam geometry. Equation (7) is originally understood as two ordered one-dimensional (1D) integrals with respect to l and θ [3, 4]. If f (x, y) and μ(x, y) are smooth, (7) is also valid in the sense of 2D integral with respect to (l, θ). With Fubini’s theorem, there would be no difference between the two ordered 1D integrals and the 2D integral, so that we are able to perform the 2D coordinate transformation and to exchange the order of the 1D integrals in (7) with respect to l and θ. The detailed derivation of (7) can be found in [3, 4]. Notice that (7) is slightly different from its original form in [3, 4]. The reason that we use the operators T and H is to emphasize the procedure of a classical FBP-type reconstruction. The attenuation map μ(x, y) is assumed to be known in this paper, thus a_θ (s, t) and h(s, θ) should be available in any coordinates because they can be calculated in advance.

Assume that the focal point positions are distributed on a circle with its center at the origin. One view of the fan-beam data acquisition geometry is described in Figure 1.

Figure 1 Each projection ray is uniquely determined by (σ, β). Points O and S are the origin and one fan-beam focal point, respectively. P is any point on which the density of radiotracer needs to be reconstructed. Here σ denotes the angle (more ...)

First, we introduce two notations of (σ, β) which will be called as fan-beam projection coordinate, and (ρ,σ)_β which denotes the polar coordinate with the origin at S. Compared with the parallel-beam geometry, (σ, β) corresponds to (l, θ) by the following relations

(8)

The Jacobian of the transformation from (σ, β) to (l, θ) is D cosσ, which is always positive for − π/2 < σ < π/2. If the support of f (x, y) and μ(x, y) is in a disk with radius less than D, the mapping between (σ, β) to (l, θ) is one-to-one inside that disk. Also, we have the following geometric relation

(9)

where

(10)

and

(11)

The detailed derivation of (9) can be found, for example, in [7, 8] by considering a simple triangular geometry. Define the weighted fan-beam projection g(σ, β) as

(12)

for the fan-beam geometry. Assume that f(x, y) and μ(x, y) are smooth, then we have

(13)

Equation (13) seems to be observed for the first time, and turns out to be the key relation to derive an FBP algorithm for fan-beam data. For each focal point position β and a fixed point (x, y), the weighting factor e^{a_θ (s, t) − h(s, θ)} is a function of the line PS. From Figure 1, the representation of PS with the coordinate of (l, θ) is (r cos(β − ), σ′ + β). Thus, in the fan-beam coordinate system (σ, β), the weighting factors in (7) can be expressed as follows,

(14)

(15)

Now, we rewrite the two terms of (7) as follows,

(16)

and

(17)

where

(18)

Apparently, both (16) and (17) follow the typical FBP procedure for the weighted fan-beam projections g (σ, β), thus an FBP reconstruction algorithm has been derived for attenuated fan-beam projections via Novikov’s inversion formula. Since both T and H are shift-invariant, they can be replaced by any other regularized kernel functions without changing the reconstruction procedure of an FBP-type. The calculations are based on the angular expression (σ, β) of projection rays, and can be extended to any other coordinate expression (η (σ), β) if η (σ) is a one-to-one mapping, e.g., equally spaced sampling along the detector.

Notice that g(σ, β) in our formula is the weighted fan-beam projection, instead of the attenuated Radon transform, and the attenuated projection has to be weighted by a complex function e^{h(s, θ)} before using our algorithm. Define the angular Hilbert transform, similar to [10], as

(19)

It has been known in [11] that the term [HRμ](l, θ) in h(s, θ) can be evaluated by the angular Hilbert transform of fan-beam projections. Actually, we have the following equation

(20)

where k(σ, β) = [Rμ](Dsinσ, σ + β) is the fan-beam projection of μ (x, y), and the angular Hilbert transform Γ is performed on the first variable of k(σ, β). A derivation of (20) can be found in Appendix A. Therefore the weighted fan-beam projection g(σ, β) can be calculated using k(σ, β) without extra interpolations between parallel-beam and fan-beam projection coordinates, i.e, we use the following formula

(21)

To summarize, the image f can be reconstructed by the following FBP-type formula

(22)

where K and σ′ are defined by (10) and (18), respectively. The attenuation map μ(x, y) and the reconstructed image are commonly expressed in the Cartesian coordinates. Also, it has been observed that the image reconstruction is not very sensitive to the noise in the attenuation map [6]. For the simplicity of implementation, we would rather evaluate a_θ (s, t) in the detector coordinate system (s, t)_θ even it could be expressed in the fan-beam projection coordinate (σ, β).

Equation (22) is an exact reconstruction formula that contains explicit differential operation. Like the Shepp-Logan filter in the classical FBP algorithm for parallel-beam projections [1], we use a regularized version of the operator T so that the explicit differential operation can be carried out by linear convolution. Hereafter, T (l) will stand for the regularized kernel function of T, and H (l) will stand for the kernel function of the Hilbert transform. More properties of H (l) and T (l) can be found in Appendix B. By using H (l) and T (l) in (7), with the relations (8)−(12), we have the following reconstruction formula

(23)

Notice that the approximation is due to the fact that T (l) is the regularized kernel for operator T.

It can be seen that the major computations come from the weighting factors W_β(r, , σ ′) and W_β¹ (r, , σ ′), which require the evaluation of attenuation-related functions a_θ(s, t) and h(s, θ). Since the attenuation map μ(x, y) is usually expressed in the Cartesian coordinates (x, y), in the evaluation of a_θ(s, t), μ_θ (s, t) was then calculated through the coordinate transformation between (x, y) and (s, t)_θ, which involved bilinear interpolations. On the other hand, in the evaluation of h(s, θ) for calculating g (σ, β), the fan-beam projection k(σ, β) of μ(x, y) was calculated through the coordinate transformation from (x, y) to (ρ,σ)_θ. Because only the real part is needed in the final reconstruction, we do not evaluate the terms that contribute to the imaginary part in (23). In summary, the attenuated fan-beam data inversion formula (23) can be realized by the following steps. (i) Pre-calculate the weighting factors e^{a_θ(s, t) − h(s, θ)} and , and multiply to the acquired fan-beam projection to obtain g (σ, β). (ii) Evaluate the two convolutions in (23). (iii) Evaluate the weighted fan-beam backprojection with the corresponding weighting functions W_β(r, , σ′) and W_β¹(r, , σ′).

(24)

where SL_δ(σ) denotes the Shepp-Logan filter and T_SL(σ)= SL_δ(σ) σ²/sin² σ. The derivation of (24) is given in Appendix B. The discrete form of T_SL (σ) is expressed as

(25)

where n ≠ 0 and T_SL (0) = 4.0 /(π δ²). The kernel Γ(σ) of the angular Hilbert transform Γ takes the following form

(26)

where σ ≠ 0 and Γ (0) = 0. Using the definition of H (l) as given in Appendix B, the function H (K sinσ) is expressed as

(27)

Let f¹ and f² be the first and the second terms in (23), respectively. Then

(28-1)

and

(28-2)

where

(29-1)

and

(29-2)

In computer implementation, parameters were discretized as dσ= δ, σ = nδ and σ ′ = n′ δ in calculating g_F¹ and g_F². In calculating f¹ and f², d β = Ω and β_k = kΩ. Notice that in the backprojection step, a linear interpolation was used for σ ′. The discrete counterparts of functions W_β and W_β¹ were A_k and B_k (their definitions are given below), respectively. Now, for a fixed point (x, y), the discretized versions of (28) and (29) can be described as the combination of the filtering step

(30-1)

and the subsequent backprojection step

(30-2)

where

and

where Δ is the sampling step of variable s in (s, t)_θ. Also notice that the calculations of A_k and B_k may need interpolations since μ(x, y) is only known on a discrete grid. Combining f¹ (x, y) and f² (x, y), image f (x, y) is reconstructed.

Figure 2 Phantoms: (a) emission distribution, (b) nonuniform attenuation map, (c) uniform attenuation map.

The geometric parameters of the fan-beam data acquisition are described as follows. The fan-focal length D = 2.0, and the subtending angle was 60º that fully covered the support region of the emission phantom and attenuation map. Both noise-free and noisy projection data were generated with 128 views evenly spaced over 360º on a circular orbit. At each view, there were 128 projection samples with equiangular fan-beam rays. The sampling rate used in this paper seems to be more realistic compared with that one in [10]. The image was reconstructed into a 128x128 image array. The reconstructed images are displayed in Figure 3 with the negative values set to zero. In all the displays, the images were linearly scaled so that their maximal values were scaled to 255. A profile drawn horizontally across the reconstructed image through the center is shown in Figure 4, where the true image values are shown. As expected, the attenuation has been compensated in the reconstructed image. For noise analysis, we define the signal-to-noise ratio (SNR) of an observed or estimated image S(m, n) as SNR(I, S) = L₂ (I) /L₂ (I − S), here I (m, n) is the true mean image. In the simulation, the noisy projection was generated using the Poisson deviate generator in [13], and the SNR of the noisy projection data was 8.07 with total counts of 641,972. The SNRs of the reconstructed images from noise-free and noisy data are 5.04 and 2.59, respectively. Numerical experiments with noisy projections were reported in [10], but no noise treatment methods were mentioned. In order to filter the noise in the reconstructed images, following the noise treatment method in [6], we used the median filter to estimate the mean of projection g(σ, β), and the 5-point symmetric Savitzky-Golay filter (see details in [13]) to smooth the filtered projection g_F¹(σ, β) in (29-1). The reconstructed image from filtered noisy projection is shown in Figure 3, and the SNR became 3.82. In this simulation, the applied filtering process does improve the SNR. The simple statistic of noisy projection data was derived to suggest optimal filters for parallel-beam projection data in [14], but the analysis might not be directly applied to the fan-beam data due to the irregularity of fan-beam geometry.

Figure 3 Reconstructed images: (a) from noise-free data, (b) from noisy data, (c) from filtered noisy data.

Figure 4 Profiles of the emission phantom and the reconstructed images.

We also conducted a simulation study when the attenuation coefficient was a constant value of 0.75 inside the largest ellipse in Figure 2. The reconstructed images are shown in Figure 5. The profiles drawn horizontally across the phantom and reconstructed images through the center are shown in Figure 6. In the simulation, the SNR of the noisy projection was 7.63 with total counts of 588,055. The SNRs of the reconstructed images from noise-free and noisy data were 4.83 and 2.38, respectively. After applying the median and Savitzky-Golay filters, the SNR of the reconstructed image became 3.60.

Figure 5 Reconstructed images: (a) from noise-free data, (b) from noisy data, (c) from filtered noisy data.

Figure 6 Profiles of the emission phantom and the reconstructed images.

In this paper, we derived a FBP algorithm for attenuated fan-beam SPECT data. The filtering is shift-invariant. Our fan-beam formula looks similar to the conventional fan-beam FBP algorithms [7, 8]. As a special case, we presented a shift-invariant FBP algorithm for fan-beam SPECT reconstruction with uniform attenuation compensation. Based on this observation, it seems that Tretiak-Metz’s inversion formula [2] might not directly accommodate an FBP algorithm for uniformly attenuated fan-beam projections with a shift-invariant filter. When the attenuation is neglected, i.e., μ 0, the second term of (23) vanishes and the first term of (23) reduces to the conventional fan-beam FBP algorithm [7, 8]. Computer simulations were provided for both non-uniform and uniform attenuation inside an ellipse. Similar to the result in [6] for parallel-beam projections, the second term of (23) was also found to have very minor contribution to the reconstruction.

In fact, there are many fan-beam data inversion formulae that can be derived through regularizing the operators T and H in different ways. For example, T (K sinσ) can be chosen as the following ramp filter

(31)

or other smoother ones like Hanning filter, and H (l) can be chosen as the similar kernels used in [4, 5]. More properties of T and H can be found in Appendix B.

A reconstruction formula, similar to [10], can be readily obtained through the following coordinate transformation

(32)

where

In this paper, we considered reconstruction algorithms for equiangular fan-beam detectors. It is straightforward to derive similar results for equally spaced fan-beam detectors (see Appendix C).

G. L. Zeng was supported in part by NIH grants CA100181 and EB00121, and Z. Liang was supported in part by NIH grant HL51466. The authors would like to express their deep thanks to reviewers for their constructive suggestions.

Let μ_β^pol (ρ,σ) denote the attenuation map in the polar coordinates (ρ, σ)_β, where the relationship between coordinates (ρ, σ)_β and (s, t) _{α + β} is shown in Figure 7.

For any σ, we have the following formula

(A1)

where the coordinate transformation from (ρ, σ)_β to (s, t)_{α + β} is

(A2)

Thus, equation (20) is derived.

Let T (l) and H (l) be the kernel functions of operators T and H, respectively. In their exact forms, T (l) and H (l) can be understood as the generalized Fourier transforms

(B1)

(B2)

It is easy to see that H₂ = −I, here I is the identity operator. For any weighting function w(l) with c₁ ≤ w(l) ≤ c₂, here c₁ and c₂ are positive constants, we have the relations

(B3)

(B4)

Let , and using the Shepp-Logan filter, we have the following regularized expression of T (K sin l) when − π/2 < l < π/2

(B5)

where is the Shepp-Logan filter, and its discrete version takes the following form,

(B6)

Therefore, equation (24) has been established.

Equation (32) is another variation of T (l) based on the use of the ramp filter. In practice, w(l) can be chosen according to the detector geometry or reconstruction factors such as the non-uniform sampling step and noise level. On the other hand, the weighting function w(l) could use different cutoff frequencies for different values l so that certain non-uniformity is included in the discretization of the exact reconstruction formula.

An imaging geometry with an equally spaced fan-beam collimator detector is shown in Figure 8, where the u-axis is parallel to the detector. The relation between the parallel-beam geometry and the fan-beam geometry is expressed as

(C1)

the Jacobian of the transformation from (u, β) to (l, θ) is D³/(D² + u²)^3/2, which is always positive.

From [12], we have the following relation

(C2)

where

(C3)

and

(C4)

Let k(u, β) represent the fan-beam projection of the same attenuation map μ(x, y) in the coordinate system (u, β), then the angular Hilbert transform Γ can be expressed as

(C5)

Therefore the weighting function, i.e., the counterpart of (20), for the equally spaced fan-beam collimation geometry can be calculated as

(C6)

Define g(u, β) as

(C7)

Following the similar calculations for obtaining equation (23), and define

(C8)

then we are able to have the following FBP reconstruction formula in the coordinate system (u, β),

(C9)

where

(C10)

and

(C11)