Search tips
Search criteria 


Logo of nihpaAbout Author manuscriptsSubmit a manuscriptHHS Public Access; Author Manuscript; Accepted for publication in peer reviewed journal;
IEEE Trans Med Imaging. Author manuscript; available in PMC 2010 May 16.
Published in final edited form as:
PMCID: PMC2871258

Reconstruction from Uniformly Attenuated SPECT Projection Data Using the DBH Method

Qiu Huang, Jiangsheng You, Gengsheng L. Zeng, Senior Member, IEEE, and Grant T. Gullberg, Fellow, IEEE


An algorithm was developed for the two-dimensional (2D) reconstruction of truncated and non-truncated uniformly attenuated data acquired from single photon emission computed tomography (SPECT). The algorithm is able to reconstruct data from half-scan (180°) and short-scan (180°+fan angle) acquisitions for parallel- and fan-beam geometries, respectively, as well as data from full-scan (360°) acquisitions. The algorithm is a DBH method, which involves the backprojection of differentiated projection data followed by an inversion of the finite weighted Hilbert transform. The kernel of the inverse weighted Hilbert transform is solved numerically using matrix inversion. Numerical simulations confirm that the DBH method provides accurate reconstructions from half-scan and short-scan data, even when there is truncation. However, as the attenuation increases, finer data sampling is required.

Keywords: Image reconstruction, attenuation, SPECT, truncation, fan-beam, half-scan, short-scan

I. Introduction

Quantitative studies in single photon emission computed tomography (SPECT) are important in clinical diagnosis and treatment. However, there are problems which can make it difficult to perform quantitative SPECT imaging. Attenuation is a major degradation factor in image quality [1]. Truncated measurements, which occur when the camera field of view (FOV) is smaller than the patient imaged, cause artifacts in the reconstructed image [2, 3]. Also, to minimize the loss in camera resolution as the source-to-detector distance increases, there is considerable interest to perform acquisitions close to the organ of interest such as in breast imaging where anterior scanning close to the breast is desired. Another advantage of anterior scanning is that a half-scan (180°) or a short-scan (180°+fan angle) can be implemented that covers the least attenuated half circles. Here we investigate a method that reconstructs the image from uniformly attenuated data measured over a half-scan acquisition for parallel-beam geometry and a short-scan acquisition for fan-beam geometry, and allows use of a relatively small detector with projections truncated at some view angles.

In the past three decades, much research has been performed in developing iterative and analytical attenuation compensation algorithms for different data acquisition geometries. Iterative algorithms, such as those in [4-12], have the potential to utilize better models of the imaging system that include various physical factors and complicated data acquisition geometries but require significant amount of computation compared to analytical methods. Whereas, analytical methods, such as those in [13-15], treat the projection data as perfect line integrals and the reconstruction can be carried out efficiently using mathematically explicit formulas. The developments of analytical methods have led to new insights into the mathematics of computerized tomography in SPECT. These insights have spurred the development of new reconstruction algorithms for various collimator geometries and detector orbit configurations with both uniform [13-15] and non-uniform [16] attenuation correction. Analytical algorithms also provide theoretical consistency conditions. These conditions are extremely important in the evaluation and understanding of the performance of iterative reconstruction algorithms and for determining if they may provide accurate solutions. Beyond the theoretical aspects there are practical reasons for developing analytical reconstruction algorithms. Analytical methods provide: a) a quick reconstruction that allows evaluation of subject positioning and determination of whether or not an experiment needs to be redone, b) a method for analytical evaluation of various factors that affect resolution and signal-to-noise, which are difficult to do with a non-linear iterative algorithm, c) a means by which to compare and to evaluate the development of iterative reconstruction methods, and d) a method useful in quality control; for example, the evaluation of whether or not geometrical calibrations, such as the center of rotation, are measured accurately.

There is considerable interest in developing algorithms that are able to reconstruct data acquired for less than 360° because of potential clinical applications. The data from a half-scan acquisition is sufficient for an exact reconstruction for the attenuation-free case [13]. However, analytical methods [14, 15] developed earlier for uniform attenuation required the projection data to be available over a full-scan. The same requirement was necessary for other analytical reconstruction methods such as the Fourier methods in [17-22], the integral geometry method in [23], and the Cormack-type inversion methods in [24, 25]. Recently progress has been made in the reconstruction of attenuated half-scan projection data [26-30]. The uniqueness of analytical reconstruction solutions when the scanning range is an open subset of [0, 360°) was proven in [31, 32], and stable reconstructions for half-scan data were later proposed for iterative algorithms in [26, 27] and an analytical algorithm in [28]. The method in [28] was an explicit formula of the convolution-backprojection type and was designed to reconstruct images from non-truncated measurements in the parallel-beam geometry. Some progress has been made in developing algorithms which are capable of dealing with fan-beam data and truncated measurements [29, 30]. The method in [29] performs least-squares fitting to find a kernel function numerically and the algorithm in [30] utilizes the calculation of a power series expansion. The algorithms in [28-30] can be implemented in a common procedure of derivative, backprojection, and Hilbert transform (DBH). In the DBH method, one of the key steps is the inversion of the finite Hilbert transform, which is an old research topic surveyed in [33]. The recent works [34-36] show progress on that subject, aiming at an exact image reconstruction from data acquisitions of less than 360°.

The DBH method has been used in x-ray CT to reconstruct a region-of-interest (ROI) from truncated data [37-39]. Recently, progress has been made [28, 30] in applying the DBH method in the reconstruction of truncated uniformly attenuated SPECT data either obtained by half-scan (parallel geometry) or short-scan (fan-beam) acquisitions. This paper further explores the DBH method for the analytical reconstruction of exponential Radon transform (eRT) data acquisitions of less than 360° in both parallel- and fan-beam acquisition geometries and develops a stable method to obtain a numerical kernel for the inverse of the finite weighted Hilbert transform. We demonstrate that the DBH method is effective in reconstructing an ROI from truncated half-scan (parallel-beam) or short-scan (fan-beam) uniformly attenuated data.

The paper is organized in five sections. In Section I, we motivate the problem under consideration. In Section II, the DBH method is developed for the 2D parallel-beam and fan-beam geometries. Section III develops the numerical implementation of the finite weighted Hilbert transform. Section IV gives implementation details and numerical results demonstrating the effectiveness of the DBH method. In Section V, we offer our conclusions and directions for further work.

II. DBH Method for Inverse of the Exponential Radon Transform

For a transaxial slice, let f(x, y) represent the distribution of a radiopharmaceutical in body tissues, which is assumed to be a smooth and compactly supported function of R2. The SPECT image reconstruction estimates f(x, y) from the detected photon counts. We denote r = (x, y) and D = {(x, y) [set membership] R2 : x2 + y2 ≤1}. We assume f(x, y) [equivalent] 0 outside of D and the attenuation of the body tissues is uniform inside the elliptical grey region in Fig. 1. For example, these assumptions can be met in brain SPECT after compensating for the attenuation of the skull [40]. We denote [theta w/ right arrow above] = (cos θ, sin θ) and [theta w/ right arrow above][perpendicular] = (−sin θ, cos θ). Let μ > 0 be the constant attenuation coefficient and Lθ,s be the distance between the point C on the boundary of the elliptical attenuator and the s-axis in Fig. 1. All photons detected in the collimator hole at (θ, s) are proportional to f(sθ+tθ)eμ(Lθ,st)dt; here the integral is carried out along the line PC. The factor e−μLθ,s can be estimated if the boundary of the uniform attenuation region is known. After multiplying the measured projections by eμLθ,s, the modified projections can be expressed as the eRT


For simplicity, we use p(θ, s) = [Rμ f](θ, s) throughout this paper. The reconstruction process estimates f(x, y) from p(θ, s).

Fig. 1
Illustration of projection coordinates corresponding to a parallel-beam scanning geometry for SPECT.

The DBH method involves calculating the derivative of the projection data, backprojection, and the inversion of the finite weighted Hilbert transform. The first two steps, the operation of derivative and backprojection, are also named as the DBP operation in [37] for the reconstruction of CT data acquired over less than 360°. According to [28], the DBP operation of the eRT is the one-dimensional (1D) weighted Hilbert transform of the original image along certain lines. Through coordinate transformations, we obtain a similar relation for the uniformly attenuated fan-beam projection data. Based on work in [30], an exact reconstruction can be achieved in a well-defined subset of the support of f(x, y) from the DBP of truncated projections.

A. DBP operation and the weighted Hilbert transform

We define an intermediate function f(x, y) as


where p′(θ, s) is the partial derivative of p(θ, s) with respect to the variable s. Equation (2) is usually called the operation of derivative and backprojection (DBP). The DBP operation yields an intermediate function f(x, y), which is different from the desired function f(x, y). Let [nabla] = ([partial differential] / [partial differential]x, [partial differential] / [partial differential]y), then from (1) we have


. Notice that by using the chain rule, we obtain


By changing the variables τ = tr · [theta w/ right arrow above][perpendicular] and using (4), we can rewrite f(x, y) as


This represents f(x, y) as the hyperbolic-cosine-weighted Hilbert transform of f(x, y) along horizontal lines. It follows that the reconstruction of f(x, y) reduces to the inversion of (5). For μ = 0, this is called the inversion of the Hilbert transform [33]. For μ ≠ 0, several formulas have recently been reported in [28, 30] to obtain f(x, y) from (5). In this paper, we refer to the DBP operation followed by the inversion of a finite weighted Hilbert transform as the DBH method for the exponential Radon transform. Note that the term “finite” refers to the limits of the last integration in (5) as being finite.

Due to the uniform attenuation, the derivative with respect to s is, in general, no longer odd, i.e., p′(θ, r · [theta w/ right arrow above]) ≠ −p′(θ + π, −r · [theta w/ right arrow above]). However, the weighted backprojection of the derivative over 360° is still zero. This property can be seen from (5) by replacing (−π / 2, π / 2), the interval of integration, with (0,2π). We also prove for a point source in Appendix A that the weighted backprojection of the derivative over 360° is zero. Due to the superposition property of a linear algorithm, the weighted backprojection for any arbitrary object is zero.

B. DBP operation for fan-beam data

A typical fan-beam data acquisition geometry with a circular focal-point trajectory is shown in Fig. 2, where each projection ray is represented by (β, σ). One particular projection ray is the arrowed line from the focal point S for the angle β with the ray angle σ. In this paper, the fan-beam uniformly attenuated projection of the function f(x, y) is defined as


where Dμf is the projection operator for the uniformly attenuated fan-beam projection data, σ [set membership] [−σm, σm], and alpha(β, σ) is a unit vector in R2 representing the direction from the focal point to the collimation hole, as shown in Fig. 2. Here, σm [set membership] (0, π / 2) denotes the maximum angle subtended by the fan-beam. Let R be the radius of the circular focal point trajectory. The coordinates in parallel- and fan-beam geometries are related by


The distance between the focal point and the s-axis is R cos σ. Thus, the definitions of (1) and (6) can be related by the following equation


We will use g(β, σ) for the modified projection of (8). By the chain rule, we have gβ(β,σ)=pθ(θ,s) and


Then, we obtain

Fig. 2
Illustration of a typical fan-beam acquisition geometry. Each projection ray is uniquely determined by (β, σ). Points O and S represent the coordinate origin and focal point, respectively. P is a point at which the value needs to be reconstructed, ...

Let r = |r| and r = (x, y) = r(cos [var phi], sin [var phi]), and K(r, [var phi], β) denotes the length of SP, and σ′(r, [var phi], β) the angle between SO and SP. It is straightforward to derive the following geometric identities [41]:




For a given r, θ, and β the following differential relationship is satisfied [42]:


For simplicity, only the fan-beam short-scan reconstruction will be discussed in this subsection. The short-scan is the scanning range of [−0.5πσm, 0.5π + σm], which is 180° plus the fan angle. We assume R sins σm = 1 so that the inner disk in Fig. 3 is the unit disk D, i.e., the projection is not truncated. In Fig. 3 for y [set membership] (−1, 1), β(y) = arcsin(y / R). After transforming the integration over θ in (2) to the integration over β in the fan-beam geometry, (2) becomes:

Fig. 3
Horizontal lines in the inner unit disk for performing the inverse of the finite Hilbert transform in the fan-beam short-scan.

The interval of integration in (15) is a subset of [−0.5πσm, 0.5π + σm], as shown in Fig. 3. Actually, for each vertical position y, the projections from β[set membership][−0.5πβ(y), 0.5π + β(y)] are sufficient to perform the DBP operation, which is similar to the π - line path integral in [43, 44]. Since the support of f(x, y) is confined inside the inner disk of Fig. 3, the DBP operation only needs to be performed on the inner disk. In other words, the function f(x, y) can be completely reconstructed in this region from uniformly attenuated fan-beam short-scan projection data.

The weight 1/ K(r,[var phi],β) may have a singularity in (15) when rR [42]. We express σ′(r,[var phi],β) in (15) in parallel-beam coordinates as


The first equality in (16) is due to (8). We modify the factor 1/ K(r,[var phi],β) in (15) using (7) and (16) to obtain


Equation (17) involves the operations of derivative and backprojection for the modified attenuated projection in a fan-beam geometry and can be readily obtained from fan-beam measurements. After (17) is implemented, the same procedure used for the parallel-beam geometry is performed to invert the finite weighted Hilbert transform.

C. Consideration of truncated data

In the previous two subsections, we assumed that the detector is large enough so that both p(θ, s) and g(β,σ) cover the support region D of f(x, y). However in SPECT imaging there are cases in which the FOV of the detector does not cover the entire object. Fig. 4 gives a typical example of such a scan where the FOV inside the thick circle does not cover the entire elliptical object. The projection acquired at the position indicated in the figure is truncated due to a small detector.

Fig. 4
Illustration of a small detector and the area of exact reconstruction by the DBH method.

Summarizing the preceding analysis, the derivative and backprojection are local operations and do not require projections to be available in the entire FOV. If f(x, y) has a support in the finite region [−L(y), L(y)], as illustrated in Fig. 4; regardless of the data acquisition geometry, f(x, y) and f(x, y) can be related by


Because of the limited data availability and the limitation in conducting the backprojection, f(x, y) may be available only in a finite region for the variable x. Intuitively, that region must be large enough to have a stable inversion. In [28], f(x, y) must be computed in (−3L(y), 3L(y)). The inversion procedure in [30] uses f(x, y) over [−L(y), L(y)] and [p(π / 2, s) + p(−π / 2, s)] / 2. The inversion algorithm presented in this work only needs f(x, y) in [−L(y), L(y)]. All those inversion formulas are mathematically exact, thus f(x, y) can be completely reconstructed from f(x, y) on [−L(y), L(y)]. Also, there is a direct extension to reconstruction of fan-beam data, as indicated in (4), (14), and (17).

III. Inversion of the Finite Weighted Hilbert Transform

It has been shown that the function f(x, y) can be obtained from the acquisition of both parallel- and fan-beam data and that f(x, y) and f(x, y) are related by a weighted Hilbert transform. In this section, we show that f(x, y) can be exactly reconstructed from f(x, y) if f(x, y) is available in the same support region of f(x, y). Frequently in SPECT imaging the projections are truncated because the FOV is not sufficiently large enough to image the entire body. However, a smaller region instead of the whole object can be accurately reconstructed from truncated projection data. In the presence of truncation, a subset of f(x, y), the region between the two dashed lines in Fig. 4, can be exactly reconstruction.

A function h(t) and its Hilbert transform H(s) are related by:


It is assumed that all singular integrals are equal to the Cauchy principal value. The function h(t) is assumed to be continuously smooth with compact support in [−q, q], q > 0. With the hyperbolic cosine weighting function, the finite weighted Hilbert transform is defined as


It is obvious that when μ = 0, the finite weighted Hilbert transform reduces to the finite Hilbert transform. Notice that Hμ (s) quickly tends to ∞ as s → ±∞.

Two classic inversion formulas [33] of the finite Hilbert transform (μ = 0) are



In order to obtain an inversion formula of the finite weighted Hilbert transform (μ ≠ 0), we construct an intermediate function ĥ(t) similar to (21), but replacing H(s) with Hμ(s) :


Upon substituting (19) into the right side of this equation and using the inversion formula in (21), we have the following series of equations that obtain a relationship between Hμ(s) and h(t)





The second equality in (23) is a result of the equality qq1[(st)q2s2]ds=0 for t[set membership](−q,q), see [33, p.174] or the alternative proof given in Appendix B. By adding the term with a value zero to Kμ(t, p), we change the singularity at s = t to a removable singularity. Notice that the function [Aμ(sp) − Aμ(tp)]/(st) is now smooth. Smoothness is important in the implementation of the algorithm. It is hard to sample functions with discontinuities, singularities, or sharp changes. Generally speaking, functions with higher frequencies need finer sampling. The smooth integrand for the definite integral makes the evaluation stable. This is desirable for efficient and accurate implementation.

It then follows that Kμ(t, p) constructs a compact operator, denoted by Ψ. Symbolically, (22) becomes


It was shown in [45] that if the eigenvalues of Ψ are different from −1, the operator (I + Ψ)−1 should exist and be bounded. We do not have a proof showing that all the eigenvalues of Ψ are different from −1 for all μ. In our simulation studies in the following section, the inverse operator (I + Ψ)−1 appears to be well defined for μ ≤ 6 when q = 1. In other words, for an organ with a 20 cm diameter, the inversion should exist for an attenuator as large as μ = 6×2 /(20 cm) = 0.6 cm−1. This is sufficiently large enough to cover the linear attenuation coefficients of most body tissues imaged in the clinic with SPECT. Hence assuming (I + Ψ)−1 exists for realistic μ values using SPECT, the original function h(t) can be completely recovered using the following inversion formula:


The expression in (26) with the dots is an operator. This operator has a functional expression when the dots are replaced with the variable t.

Equation (25) is a Fredholm integral equation of the second kind. According to the comments in [46], the matrix inversion based algorithm discussed in [47, 48] is a stable way to numerically solve (25). Using the trapezoidal rule for evenly sampled data, the numerical implementation of (26) is summarized in the following four steps.

Step1: Initialize the following

  • Compact support interval is [−1, 1], i.e., q = 1.
  • Sampling interval is Δt = 1/M with M = 512 for two 1D signals and M = 128 for SPECT data in this paper.
  • Sampling points within the support are tm = (m + 0.5)Δt and −MmM −1.

Step 2: Compute the kernel Kμ(t, p)

  • Aμ(t)={(cosh(μt)1)t,|t|0.5Δt0,|t|<0.5Δt}.
  • Mμ(s,t,p)={Aμ(sp)Aμ(tp)π2(st),|st|0.5Δt0.5(μπ),|st|<0.5Δt}.
  • Kμ(tm,tn)=Δt1tm2Σj=MM1[Mμ(tj,tm,tn)1tj2].
  • (I + Ψ)−1 = {Kμ(tm, tn) + δm,n}−1, where δm,n is the Kronecker delta function.

Step 3: Perform the integral transform


In the implementation, the contribution of the integrand at singular points is not sampled.

Step 4: Invert the matrix (I + Ψ)


The calculation of Kμ(tm, tn) in step 2 does not involve any singularity; thus, all numerical computations are stable. When applying this numerical procedure to the image reconstruction problem (5), we let Hμ(s) = f(x = s, y) and find f(x, y) = h(t = x) for a fixed y. This row-by-row finite inversion procedure is repeated for all y coordinates.

IV. Computer Simulation Results

The key step in the DBH method is the inversion of the finite weighted Hilbert transform. The numerical behavior of (19) and (26) was first studied, and then the DBH method was applied to the inversion of the eRT for half-scan parallel-beam and short-scan fan-beam data acquisition geometries both with truncated and non-truncated projections.

A. Numerical realization of (19) and (26) for 1D signals

Two 1D signals h(t) = sin(πt) and h(t) = 1− | t | with | t |< 1 were chosen to investigate the numerical behaviors of the finite weighted Hilbert transform (19) and the inversion formula (26). Original signals were evenly sampled over 1024 points for | t |< 1 as described previously in the implementation steps in Section III. A very high sampling rate was used in this simulation to ensure that the finite weighted Hilbert transform Hμ(s) was accurate. For this simulation, the attenuation coefficient μ = 6 was equivalent to 6/512 per sampling interval. If 1024 points are sampled along an object diameter of 20 cm, this is equivalent to an attenuation coefficient equal to 0.6 cm−1. In this simulation, the weighted Hilbert transform of the original signal was calculated using (19) and is shown in Figs. 5(a) and 5(b). Then these two 1D signals were reconstructed from their finite weighted Hilbert transform according to the four steps described in Section III. For comparison, the original signals and their reconstructions are plotted in Figs. 5(c) and 5(d). Numerically, equation (26) does provide fairly accurate reconstructions.

Fig. 5
Weighted Hilbert transform Hμ(s) for (a) h(t) = sin(πt) and (b) h(t) = 1− | t |, only shown for |s|<1. Original signals and their reconstruction from the finite weighted Hilbert transform: (c) for h(t) = sin(πt ...

B. Parallel-beam half-scan with and without truncation and with and without noise

In the SPECT simulation study, the object function f(x, y) is chosen to be the modified Shepp-Logan phantom shown in Fig. 6(a), which is composed of 10 ellipses. The intensity in each ellipse has been modified from the original Shepp-Logan phantom [13] to mimic the distribution of a radiopharmaceutical concentration. These ellipses were described in [30], and were scaled to be contained in the unit disk D in our simulations. The phantom is represented by a 256×256 grid. In terms of pixels, the lengths of the major and minor axes of the most outer ellipse are 184 and 138, respectively. The three detectors around the phantom indicate the range of the projections being sampled. Projection data were evenly sampled on a grid of 400×256 over [π / 2, 3π / 2)×(−1, 1), and are shown in Fig. 6(b). All simulated projection data were analytically generated based on the weighted line integrals of the modified Shepp-Logan phantom with a uniform attenuation coefficient of 0.15 cm−1.

Fig. 6
(a) The modified Shepp-Logan phantom. The numbers in the figure without dimensions indicate the relative distributions of the radiopharmaceutical for each of the ellipses. For example, the five small ellipses have a value 0.4. These values were scaled ...

In the study of the effect of noise with the DBH method, we converted the analytically generated projection data into Poisson random samples. The total counts were 2×107, with each bin containing about 200 photons on average. Projection data were first modified to place the detector at the center of rotation, as explained in Section II. Once the modified projections p(θ, s) were obtained, the DBP image f(x, y) was calculated by (2) on a 256×256 grid, the same size used to sample the original modified Shepp-Logan phantom. At each fixed view angle, the modified projections were convolved with a discrete band-limited five-point kernel to yield the derivative p′(θ, s). The numerical implementation of the derivative was a mid-point method. The DBP image was obtained by backprojecting the derivative over 180° with an exponential weighting function eμr·θ[perpendicular]. Then the finite inverse Hilbert transform (26) was performed to find f(x, y) along each row. In other words, using the four steps in Section III, the desired image f(x, y) (corresponding to h(t = x)) at each fixed y is reconstructed from its weighted Hilbert transform f(x, y) (corresponding to Hμ (s) = f(x = s, y)) at the same y.

The reconstructed images from noise-free and noisy projection data are shown in Fig. 7(a). The grayscale is [0.15, 0.45] for the display of all the reconstructed images throughout the paper. Two profiles passing through the reconstructed images at y = 0 are plotted in Fig. 7(b) to compare the difference between the reconstruction and the original phantom. For the noise-free data, the visual observation from Fig. 7(a) and the 1D profiles in Fig. 7(b) show that the reconstruction is fairly accurate.

Fig. 7
(a) Reconstructed images from parallel-beam half-scan, noise-free (left), and noisy (right) data (total counts equal 2×107). The grayscale is [0.15, 0.45]. (b) Profiles of reconstructed images from half-scan, noise-free, and noisy data (total ...

For the noisy data, there is a non-uniform distribution of the noise from the top to the bottom in the reconstructed image. More specifically, the upper half of the reconstructed image looks noisier than the lower half. For the case when the total counts equal 2×107, the percent root mean square errors for the two circles shown in the right image of Fig. 7(a) are 7.33% for the lower and 7.67% for the upper circle. Here the root mean square error was defined as the ratio of the standard deviation and the mean value sampled over the corresponding regions of interest.

The difference in the noise structure between the upper and lower half of the image is related to the weighting factor eμr · [theta w/ right arrow above][perpendicular] in the backprojection. If r · [theta w/ right arrow above][perpendicular] > 0, the weighting factor eμr · [theta w/ right arrow above][perpendicular] is less than 1 and suppresses the noise in the backprojection; otherwise, it magnifies the projection noise. For our simulations, the projections were sampled from θ [set membership][π / 2, 3π / 2]. Therefore, the backprojection of the differentiated projections are weighted by a factor greater than 1 in the upper half of the image because r · [theta w/ right arrow above][perpendicular] < 0 and less than 1 in the lower half of the image because r · [theta w/ right arrow above][perpendicular] > 0. Thus, the reconstruction in the upper half of the image is noisier.

For those pixels being generally further from the detector, this effect of the noise being worse can be observed from the variance images for three different noise levels in Fig. 8. For total counts of 1×105, 5×105 and 1×106, we performed 1000 simulations without truncated projections with random noise generated for each realization. The ensemble variances [49, 50] of the reconstructed images were calculated and are shown in Fig. 8 for the three noise levels. The brightness at each pixel is the value of the variance at that pixel. The images in Fig. 8 show higher variance in the top-half of each image. Also, it is clear that the higher the total counts, the lower the ensemble variance. The mean and the standard deviation of the value for the pixel at the origin (x, y) = (0,0) are tabulated in Table 1. The Tests 1, 2, and 3 refer to the tests performed for total counts of 1×105, 5×105 and 1×106, respectively. When the counts are low, we observe that the noise causes large deviations. In order to investigate how well the sample mean estimates the population mean we calculated the 95% confidence interval (Xσnt,X+σnt) for the sample mean X with sample standard deviation σ [51]. For n = 1000 realizations, the degrees of freedom is 999 giving a critical value of t = t0.025,n−1 = 1.962 for the t-distribution from the table on page 485 of [51]. The lower and upper confidence limits were calculated and are tabulated in the last two columns of Table 1. As can be seen in Fig. 9, the confidence interval shrinks as the noise decreases, which means the estimation is more precise.

Fig. 8
Variance images for different noise levels without truncation. The total counts equal 1×105, 5×105 and 1×106, from left to right. The grayscale is [0, 0.1155] for all the three images.
Fig. 9
The 95% confidence interval for the pixel value at the origin. The horizontal axis indicates six tests: Tests 1, 2, 3 for non-truncated data with total counts of 1×105, 5×105 and 1×106, and Tests 4, 5, 6 for truncated data with ...
Table 1
Statistics used in the calculation of the 95% confidence interval for the sample means. Tests 1, 2, 3 are for non-truncated data with total counts of 1×105, 5×105 and 1×106, and Tests 4, 5, 6 are for truncated data with total counts ...

For the study with truncated data, we considered the detector to sample only the central 156 rays of the modified Shepp-Logan phantom at each angle as illustrated for the elliptical object in Fig. 4. This small detector defines the FOV to be within the thick circle with a radius of 78. Recall that the maximum length of the ellipse of the modified Shepp-Logan phantom was 184; thus, a total of 28 pixels along the major axis were truncated. As pointed out previously, the derivative and backprojection are local operations; thus, the values of the intermediate function f(x, y) between the two dashed lines illustrated in Fig. 4 are sufficient to provide an exact reconstruction of f(x, y) within that region. In this simulation, 400×156 noise-free and noisy truncated projection data were generated, and (2) was then used to backproject the derivative of the truncated data. The inverse (26) of the finite weighted Hilbert transform was applied to reconstruct the object function f(x, y) between the two dashed lines in Fig. 4.

The reconstructed images are shown in Fig. 10. Three variance images (see Fig. 11) were calculated from 1000 noise realizations with truncation for total counts of 1×105, 5×105 and 1×106, respectively. The 95% confidence intervals for the mean value of the pixel at the origin (x, y) = (0,0) are tabulated in Table 1. The Tests 4, 5, and 6 refer to the tests performed for the truncated case for total counts of 1×105, 5×105 and 1×106, respectively. The same trend as in the non-truncated case can be observed. The higher the total number of counts, the lower the ensemble variance and the narrower the confidence interval. In other words, greater precision is expected when we increase the counts. This is shown graphically in Fig. 9. Note that the true value at the origin was 0.3. The results show differences between the reconstruction of truncated and non-truncated data in both the bias and the variance. For the non-truncated cases the bias and precision improve with the increase in the total counts whereas for the truncation cases the precision improves with increased counts but the bias appears to become worse.

Fig. 10
Reconstructed images from parallel-beam, truncated noise-free (left), and noisy (right) data. The grayscale is [0.15, 0.45].
Fig. 11
Variance images for different noise levels with truncation. The total counts equal 1×105, 5×105 and 1×106, from left to right. The grayscale is [0, 0.0339] for all the three images.

C. Fan-beam truncated short-scan data with and without noise

Equation (26) indicates that the data range required for Hμ (s) is the same interval as for h(t), indicating that the DBP image f(x, y) in the inner disk of Fig. 3 is sufficient for reconstructing the original image f(x, y). In order to ensure a numerically stable inversion using (26), the DBH operation has to be conducted in a larger area than the support of f(x, y). Since the fan-beam short-scan reconstruction using the DBH method with and without truncation is virtually identical within the truncated FOV, we only show in this paper the simulation for the reconstruction of truncated short-scan data. In our computer simulations, the focal length was R = 400 pixels with a subtending angle of σm = 32.5°. The modified data g(β,σ) were analytically computed on a 256×256 grid for [−122.5°, 122.5°]×[−32.5°, 32.5°]. We also converted g(β,σ) into Poisson random samples to simulate the noisy projections for total counts equal to 2×107.

The derivative of a sampled projection was calculated using (10), the backprojection was performed by strictly following (15) inside the thick rectangle illustrated in Fig. 12, and then (26) was applied to reconstruct f(x, y). The reconstructed images from noise-free and noisy data are shown in Fig. 13. The direction for the finite inverse weighted Hilbert transform was set to be horizontal as described in Fig. 3. To reconstruct the value at point a, data over the arc AMB were used. While at point b, only data over the arc CMD were used. The different ranges of the backprojection in (15) over β for different points such as point a and point b cause different sampling densities. Such a nonstationary sampling density can be mitigated through the parallel-beam coordinate expression (17), which was suggested in [42] for the reconstruction of fan-beam CT data that is partially sampled.

Fig. 12
The outer circle is the trajectory of the focal point for the fan-beam geometry, and the inner circle is the support region of the object function f(x, y). The DBP operation is performed inside the rectangular area.
Fig. 13
Reconstructions of truncated fan-beam, noise-free (left), and noisy (right) data. The grayscale is [0.15, 0.45].

In subsection A of section IV, a very fine sampling was used to implement both (19) and (26) for two reasons: one was to calculate Hμ(s) with high accuracy and the other was to reduce the numerical errors in the computation of (I + Ψ)−1. In subsections B and C, a sampling of 400 projection angles over π with each projection having 256 samples, was used for sampling p(θ, s) or g(β,σ) in order to make sure the weighted Hilbert transformed image f(x, y) was accurate enough. In our numerical simulation study, it was found that the sampling is related to the attenuation coefficient μ with finer sampling required for higher attenuation coefficients. In subsection C of section IV, we calculated the derivative of the projection data using (10). However, the derivative with respect to β can be avoided by using the following equation:


Using (27) could probably help reduce the numerical errors in calculating the partial derivative with respect to β, which is usually under-sampled in SPECT imaging.

V. Discussion and Conclusion

The DBH method was applied to the inversion of the eRT for both parallel- and fan-beam data acquisition geometries for SPECT. In the DBH method, the backprojection of the differentiated projections was obtained first and an inverse of the finite weighted Hilbert transform was performed to restore the image. Computer simulations show that the DBH algorithm produces a reconstruction of uniformly attenuated SPECT data when the attenuation coefficient is 0.15cm−1. The effectiveness of our method was demonstrated with the reconstruction of parallel-beam half-scan and fan-beam short-scan truncated projections. However, issues such as sampling resolution, noise effects, and clinical efficacy have to be further investigated.

We developed a matrix inversion method to numerically generate the kernel for the inverse of the finite weighted Hilbert transform. Our numerical method consists of forming a small nonsingular matrix (256×256 in our computer simulations) by numerically evaluating definite integrals in (23), which do not contain any essential singularities. The singularity at s = t in (23) was transformed into a removable singularity by adding an integral term equal to zero. This matrix does not have any zero eigenvalues. The eigenvalues depend on the attenuation coefficient. For attenuation coefficients in the range of most clinical nuclear medicine applications, our computer simulations indicate that the inversion of the matrix is stable and insensitive to computer numerical errors.

The most important aspect of this work is the development of a stable method to obtain a numerical kernel of the inverse of the finite weighted Hilbert transform. A new result in this paper is the application of the DBH algorithm with this new kernel to parallel-beam half-scan, fan-beam short-scan, and truncated data. This method requires Hμ(s) to be known in the interval (−q, q), whereas in [28] Hμ (s) is required in the interval (−3q, 3q). In [30] a power series expansion method was used to solve for the inverse to the integral equation in (19) starting with (20). Our proposed matrix inversion method started instead with (21) to derive an inversion procedure. In [29] a different matrix inversion method was developed, where the equation in (5) of this paper was first discretized and transformed into a system of linear equations, a matrix inversion method was used to solve this system. Computer simulations showed that the method developed in this paper is more stable than the one in [29].

The analytical reconstruction algorithm developed in this paper for the correction of uniform attenuation is computationally more efficient than an iterative algorithm and may have potential applications in clinical SPECT. A possible application is in brain imaging after the data have been compensated for the attenuation of the bone. The two dimensional (2D) algorithm could also be extended to the development of a fully three dimensional (3D) analytical reconstruction algorithm with uniform attenuation correction for breast imaging with rotating slant hole collimation [52]. This would be an important extension of the 2D algorithm to fully 3D tomographic reconstruction. It is also possible that the DBH reconstruction method could be applied to other more complicated geometries such as variable focal-length fan-beam geometries.

The DBH method is a new direction in the development of analytical reconstruction algorithms for tomographic imaging. It provides important insight into understanding the limitations and possibilities for reconstructing truncated projections with uniform attenuation correction. The development of this method is actively being carried out in CT [37-39] and has potential for SPECT [28-30] with uniform attenuation. There is potential of extending our algorithm to fully 3D applications and to the reconstruction of cone-beam data with correction for uniform attenuation [53]. A more difficult problem, which is still unsolved, is the development of an analytical algorithm for reconstructing an ROI from truncated projections with a non uniform attenuation distribution.


This work was supported in part by the Director, Office of Science, Office of Biological and Environmental Research, Medical Science Division of the U. S. Department of Energy under Contract No. DE-AC03-76SF00098 and in part by Public Heath Service grant R01EB000121 awarded by the National Institute of Biomedical Imaging and Bioengineering, Department of Health and Human Services.


A. Proof that the DBP of a point source over 360° is equal to zero

Without loss of generality, let the object be a delta function at (0, 1), i.e., f(x, y) = δ(x, y−1). According to the definition in (2), the eRT is p(θ, s) = eμcosθ δ(s − sin θ). The derivative of the eRT with respect to s can be expressed as p′(θ, s) = eμcosθ δ′(s − sinθ).

Next, consider the weighted backprojection of p′(θ, s) = eμcosθ δ′(s − sinθ) at a fixed point r = (x, y):


Let θ0 be the particular detector angle where the line connecting r and the point source location (0,1) is orthogonal to the detector as shown in Fig. 14.

Fig. 14
Acquisition of a point source activity.

Denote s0 = sin θ0 = r · θ0 and t0 = r · θ0[perpendicular]. Using the properties of the delta function, the above weighted backprojection can be written as


where g(θ) = r · θ − sin θ = x cos θ + y sin θ − sin θ and g′(θ) = −x sin θ + y cos θ − cos θ = r · θ[perpendicular] − cos θ.

B. An alternative proof of qq1(st)q2s2ds=0 for t [set membership] (−q, q)

Let x = st, then


This integral is evaluated in the sense of the Cauchy Principle Value. That is,


In the following we use the formula (2.266) in [54]:


For negative x let x = −y, then we have


Using the short-hand notation F(x), the integral I can be evaluated as


Since −q < t < q, (qt) > 0 and (−qt) < 0,



Therefore, F(qt) − F(−qt) = 0. Also,


Evaluating the limε→0 (F(−ε) − F(ε)), we have


Thus, I = 0.

Contributor Information

Qiu Huang, E. O. Lawrence Berkeley National Laboratory, One Cyclotron Road, Mail Stop 55R0121, Berkeley, CA 94720, USA.

Jiangsheng You, Cubic Imaging LLC, 18 Windemere Dr, Andover MA 01810, USA.

Gengsheng L. Zeng, Utah Center for Advanced Imaging Research (UCAIR), Department of Radiology, University of Utah, 729 Arapeen Drive, Salt Lake City, UT 84108, USA.

Grant T. Gullberg, E. O. Lawrence Berkeley National Laboratory, One Cyclotron Road, Mail Stop 55R0121, Berkeley, CA 94720, USA.


1. Tsui BMW, Frey EC, LaCroix KJ, Lalush DS, McCartney WH, King MA, Gullberg GT. Quantitative myocardial perfusion SPECT. J Nucl Cardiol. 1998;5:507–522. [PubMed]
2. Chen J, Galt JR, Case JA, Ye J, Cullom SJ, Durbin MK, Shao L, Garcia EV. Transmission scan truncation with small–field-of-view dedicated cardiac SPECT systems: impact and automated quality control. J Nucl Cardiol. 2005;12:567–573. [PubMed]
3. Celler A, Dixon KL, Chang Z, Blinder S, Power J, Harrop R. Problems created in attenuation-corrected SPECT images by artifacts in attenuation maps: A simulation study. J Nucl Med. 2005;46:335–343. [PubMed]
4. Shepp LA, Vardi Y. Maximum likelihood reconstruction for emission tomography. IEEE Trans. Med. Imag. 1982;Mi-1:113–122. [PubMed]
5. Lange K, Carson R. EM reconstruction algorithms for emission and transmission tomography. J. Computer Assisted Tomography. 1984;8:306–316. [PubMed]
6. Hudson HM, Larkin RS. Accelerated image reconstruction using ordered subsets of projection data. IEEE Trans. Med. Imag. 1994;13:601–609. [PubMed]
7. Gullberg GT, Huesman RH, Malko JA, Pelc NJ, Budinger TF. An attenuated projector-backprojector for iterative SPECT reconstruction. Phys Med. Biol. 1985;30:799–816. [PubMed]
8. Manglos SH, Jaszczak RJ, Floyd CE, Hahn LJ, Greer KL, Coleman RE. Non-isotropic attenuation in SPECT: Phantom test of quantitative effects and compensation techniques. J. Nucl. Med. 1987;28:1584–1591. [PubMed]
9. Tsui BMW, Hu H-B, Gilland DR, Gullberg GT. Implementation of simultaneous attenuation and detector response correction in SPECT. IEEE Trans. Nucl. Sci. 1988;35:778–783.
10. Tsui BMW, Gullberg GT, Edgerton ER, Ballard JG, Perry JR, McCartney WH, Berg J. Correction of nonuniform attenuation in cardiac SPECT imaging. J. Nucl. Med. 1989;30:497–507. [PubMed]
11. Zeng GL, Gullberg GT, Tsui BMW, Terry JA. Three dimensional iterative reconstruction algorithms with attenuation and geometric point response correction. IEEE Trans. Nucl. Sci. 1991;38:693–702.
12. Liang Z, Turkington TG, Gilland DR, Jaszczak RJ, Coleman RE. Simultaneous compensation for attenuation, scatter and detector response for SPECT reconstruction in three dimensions. Phys. Med. Biol. 1992;37:587–603. [PubMed]
13. Shepp L, Logan B. The Fourier reconstruction of a head section. IEEE Trans. Nucl. Sci. 1974;21:21–43.
14. Tretiak O, Metz C. The exponential Radon transform. SIAM J. Appl. Math. 1980;39:341–354.
15. Gullberg GT, Budinger TF. The use of filtering methods to compensate for constant attenuation in single-photon emission computed tomography. IEEE Trans. BioMed. Eng. 1981;28:142–157. [PubMed]
16. Novikov RG. An inversion formula for the attenuated X-ray transformation. Ark. Math. 2002;40:145–167.
17. Bellini S, Piacenti M, Caffario C, Rocca F. Compensation of tissue absorption in emission tomography. IEEE Trans. Acoust. Speech Signal Processing. 1979;27:213–218.
18. Markoe A. Fourier inversion of the attenuated x-ray transform. SIAM J. Math. Anal. 1984;15:718–722.
19. Inouye T, Kose K, Hasegawa A. Image reconstruction algorithm for single-photon-emission computed tomography with uniform attenuation. Phys. Med. Biol. 1989;34:299–304. [PubMed]
20. Liang Z, Ye J, Hanington DP. An analytical approach to quantitative reconstruction of nonuniform attenuated brain SPECT. Phys. Med. Biol. 1994;39:2023–2041. [PubMed]
21. Metz CE, Pan X. A unified analysis of exact methods of inverting the 2-D exponential Radon transform, with implications for noise control in SPECT. IEEE Trans. Med. Imag. 1995;17:643–658. [PubMed]
22. Liang Z, Ye J, Cheng J, Harrington DP. The inversion of the exponential Radon transform for quantitative brain SPECT. Phys. Med. Biol. 1996;41:1227–1232. [PubMed]
23. Palamodov VP. An inversion method for an attenuated x-ray transform. Inverse Problems. 1996;12:717–729.
24. Puro A. Cormack-type inversion of exponential Radon transform. Inverse Problems. 2001;17:179–188.
25. You J, Liang Z, Zeng GL. A unified reconstruction framework for both parallel-beam and variable focal-length fan-beam collimators by a Cormack-type inversion of exponential Radon transform. IEEE Trans. Med. Imag. 1999;18:59–65. [PubMed]
26. Noo F, Wagner JM. Image reconstruction in 2D SPECT with 180° acquisition. Inverse Problems. 2001;17:1357–1371.
27. Pan X, Kao CM, Metz CE. A family of π-scheme exponential Radon transforms and the uniqueness of their inverses. Inverse Problems. 2002;18:825–836.
28. Rullgård H. An explicit inversion formula for the exponential Radon transform using data from 180 degrees. Ark. Math. 2004;42:353–362.
29. Huang Q, Zeng GL, Gullberg GT. An analytical inversion of the 180° exponential Radon transform with a numerically generated kernel. International Journal of Image and Graphics. 2007;7:71–85.
30. Noo F, Defrise M, Pack J, Clackdoyle R. Image reconstruction from truncated data in SPECT with uniform attenuation. Inverse Problems. 2007;23:645–667.
31. Hazou IA, Solmon DC. Inversion of the exponential x-ray transform I: analysis. Math. Methods Appl. Sci. 1988;10:561–574.
32. Hazou IA, Solmon DC. Filtered-backprojection and the exponential Radon transform. J. Math. Anal. Appl. 1989;141:109–119.
33. Tricomi FG. Integral Equations. Interscience Publisher Inc; New York: 1957.
34. You J, Zeng GL. Explicit finite inverse Hilbert transforms. Inverse Problems. 2006;22:L7–L10.
35. Defrise M, Noo F, Clackdoyle R, Kudo H. Truncated Hilbert transform and image reconstruction from limited tomographic data. Inverse Problems. 2006;22:1037–1055.
36. Zeng GL, You J, Huang Q, Gullberg GT. Two finite inverse Hilbert transform formulae for region-of-interest tomography. International Journal of Imaging and Technology. 2007;17(Issue 4):219–223.
37. Noo F, Clackdoyle R, Pack JD. A two-step Hilbert transform method for 2D image reconstruction. Phys. Med. Bio. 2004;47:3903–3923. [PubMed]
38. Zou Y, Pan X. Exact image reconstruction on PI lines from minimum data in helical cone-beam CT. Phys. Med. Biol. 2004;49:941–959. [PubMed]
39. Clackdoyle R, Frederic N, Guo J, Roberts JA. Quantitative reconstruction from truncated projections in classical tomography. IEEE Trans. Nuc. Sci. 2004;51:2570–2578.
40. Liang Z, Ye J, Harrington DP. Three-Dimensional Image Reconstruction in Radiology and Nuclear Medicine. Kluwer Academic Publishers; 1996. Quantitative brain SPECT in three dimensions: an analytical approach without transmission scans; pp. 117–132.
41. Kak AC, Slaney M. Principles of Computerized Tomography. IEEE; Piscataway, NJ: 1987.
42. You J, Zeng GL. Hilbert transform based FBP algorithm for fan-beam CT full and partial scans. IEEE Trans. Med. Imag. 2007;26:190–199. [PubMed]
43. Katsevich AI. Analysis of an exact inversion algorithm for spiral cone-beam CT. Med. Phys. 2001;47:2583–2597. [PubMed]
44. Katsevich AI. An improved exact filtered backprojection algorithm for spiral computed tomography. Adv. Appl. Math. 2004;32:681–697.
45. Rudin W. Functional Analysis. McGraw-Hill Book Company; New York: 1973.
46. Press WH, Teukolsy SA, Vetterling WT, Flannery BP. Numerical Recipes in C. Cambridge University Press; 1992.
47. Atkinson KE. A Survey of Numerical Methods for the Solution of Fredholm Integral Equations of the Second Kind. SIAM; Philadelphia: 1976.
48. Delves LM, Mohamed JL. Computational Methods for Integral Equations. Cambridge University Press; Cambridge, UK: 1985.
49. Wilson DW, Tsui BMW. Noise properties of filtered-backprojection and ML-EM reconstructed emission tomographic images. IEEE Trans. Nucl. Sci. 1993;40:1198–1203.
50. Barrett HH. Objective assessment of image quality: effects of quantum noise and object variability. J. Opt. Soc. Am. A. 1990;7:1266–1278. [PubMed]
51. Zar JH. Biostatistical Analysis. 2nd ed Prentice Hall, New Jersey: 1984.
52. Xu J, Tsui BM. An accurate region-of-interest reconstruction method using a two-step Hilbert transform for rotating slant hole SPECT. the 54th SNM annual meeting; San Diego. June 3-6, 2006.
53. Noo F, Pack JD. Theory for image reconstruction from divergent-beam projections in SPECT; IEEE Nuclear Science Symposium Conference Record; 2006. pp. 3449–3452.
54. Gradshteyn IS, Ryzhik IM. Table of Integrals, Series, and Products. 5th Ed Academic Press; San Diego: 1994.