Abstract
We improve data extrapolation for truncated computed tomography (CT) projections by using Helgason-Ludwig (HL) consistency conditions that mathematically describe the overlap of information between projections. First, we theoretically derive a 2D Fourier representation of the HL consistency conditions from their original formulation (projection moment theorem), for both parallel-beam and fan-beam imaging geometry. The derivation result indicates that there is a zero energy region forming a double-wedge shape in 2D Fourier domain. This observation is also referred to as the Fourier property of a sinogram in the previous literature. The major benefit of this representation is that the consistency conditions can be efficiently evaluated via 2D fast Fourier transform (FFT). Then, we suggest a method that extrapolates the truncated projections with data from a uniform ellipse of which the parameters are determined by optimizing these consistency conditions. The forward projection of the optimized ellipse can be used to complete the truncation data. The proposed algorithm is evaluated using simulated data and reprojections of clinical data. Results show that the root mean square error (RMSE) is reduced substantially, compared to a state-of-the-art extrapolation method.
1. Introduction
It is known that traditional computed tomography (CT) reconstruction algorithms, for example, filtered backprojection methods, are not compatible to laterally truncated projection data, which appears often in the case of either when the object extends outside of the field of view (FOV) or X-ray beam collimation for the purpose of dose reduction. If data truncation is not effectively compensated for, it will result in cupping-like artifacts and incorrect gray-value levels in the reconstruction. A typical approach to reduce truncation artifacts is to perform extrapolation, for example, with the symmetric mirroring method [1], water cylinder extrapolation method [2], optimization-based extrapolation scheme [3], or implicit extrapolation method performed in the second-order derivative domain [4]. However, these heuristic extrapolation methods typically rely on techniques that complete the truncated data by means of a continuity assumption and thus appear to be ad hoc.
It has been demonstrated that any physically consistent sinogram has a strong restriction in its functional form [5]. This restriction is expressed by Helgason-Ludwig (HL) consistency conditions [6, 7], which are a mathematical expression to precisely describe the overlap of information between different projections. The HL consistency conditions play an important role in image reconstruction from imperfect projection data (e.g., due to noise, motion, and truncation) since these projections no longer satisfy the HL conditions. Related work uses HL conditions to estimate motion parameters directly from sinograms [8–10] or to solve the problem of limited angle tomography using a variational formulation that incorporates HL conditions [11]. In PET/SPECT, the HL consistency conditions were also used for attenuation correction if no transmission data is available [12].
This work addresses consistency-based sinogram completion. The methods proposed in [2, 13] implicitly used the zeroth-order HL consistency condition; that is, the DC term is the same for all projections, as a constraint for data extrapolation. The first-order condition, which corresponds to the first moment of the projections and describes the so-called “center of mass,” was also used to guide the extrapolation procedure [14]. Later, the elliptical extrapolation suggested in [15] explicitly used a small subset of the consistency conditions in the original HL formulation (projection moment theorem) so that large numerical instability can be avoided when computing the moment terms. The approach in [16] modified this original formulation by expanding the Radon transform in terms of its basis functions and incorporated not only one or two HL consistency conditions, but theoretically an infinite number of such constraints. However, the HL consistency conditions proposed in [16] were represented in the Chebyshev-Fourier domain, which increased computational complexity for practical applications. To simplify the computation, the method in [17] refined the Chebyshev-Fourier representation of HL conditions using an FFT with additional cosine transform along the detector channel. Furthermore, fan-beam to parallel-beam rebinning is required since the consistency conditions were only derived for parallel-beam geometry.
In this paper, we first derive the HL consistency conditions in the 2D Fourier domain from their original formulation. The Fourier representation shows that there is a zero energy region appearing in the Fourier transform of a sinogram (symmetric for parallel-beam and asymmetric for fan-beam geometry). This property was also demonstrated in [18, 19], which is referred to as the Fourier property of a sinogram and which was approximately arrived at using the parallel-/fan-beam sinogram of a delta point object. If the projection data is imperfect or incomplete, the zero energy double-wedge region contains nonnegligible values that indicate the corresponding inconsistent components. Several applications using this Fourier property of a sinogram can be found in [8, 20–23]. In this work, we theoretically prove the equivalence between the HL consistency conditions and the Fourier property of a sinogram and show the advantages of applying these Fourier-based consistency conditions: first, an infinite number of conditions are considered; and second, 2D Fourier transform via FFT is computationally more efficient than the Chebyshev-Fourier transform [16] or Lagrange-Fourier transform [11]. These features allow us to develop an efficient data extrapolation method by optimization of a cost function based on the Fourier-based HL conditions. First investigation on the method was also reported in [24].
The organization of the paper is as follows. In Section 2, we review the HL consistency conditions in its original formulation and the modified Chebyshev-Fourier representation. Then, we derive 2D Fourier-based HL consistency conditions and extend the conditions from parallel-beam to fan-beam geometry for centered objects. In Section 3, we design a cost function based on HL consistency conditions, which we use in a constrained optimization over a uniform ellipse that describes the object outline. In Section 4, we present experimental results from both a simulated phantom and reprojections of clinical data. In Sections 5 and 6, we discuss the relevant issues and draw conclusion.
2. Consistency Conditions
2.1. Helgason-Ludwig (HL) Consistency Conditions
In this section, we review the original formulation of HL consistency conditions, which is also referred to as the projection moment theorem in the literature [16]. Suppose the object is supported on the unit disk centered at the origin. Let be the th moment of the sinogram with respect to the detector bin , which is defined as
Then, the function does not change arbitrarily when the rotation angle varies. The Fourier series expansion of can be written as follows: with Fourier coefficients given by Then, it is readily proven [25] that all necessarily satisfy
2.2. Chebyshev-Fourier Representation of HL Conditions
The derivation of the Chebyshev-Fourier version of HL conditions is similar to the one in [16]. But here we use the Chebyshev polynomial of the first kind to replace the monomial term , instead of the Chebyshev polynomial of the second kind as shown in [16].
Note that the functions do not form a set of orthogonal basis functions on (the Radon transform maps the Hilbert space consisting of finite norm objects to the Hilbert space consisting of finite norm sinograms ), where . In the following we show that the monomial can be replaced by the th-order orthogonal polynomial, such that the HL consistency conditions are tractable to use in reconstruction.
The sinogram can be expanded in a series as follows:where denote the expansion coefficients and denotes the th-order Chebyshev polynomial of the first kind, which is defined bywith
Let ; we define the inner product of and as follows: From Appendix A we will prove that form an orthogonal basis of .
Then, we can obtain an expression of the expansion coefficients as a scalar product
By comparing (5) and (10), it is noted that the coefficients are related to by the following combination: Then, we will have the HL conditions as
In sum,
2.3. 2D Fourier-Based HL Consistency Conditions
We perform the 2D Fourier transform to both sides of (6)
Because the term is uniformly convergent, the order of the integral operator and the sum operator can be changed:
Because of the orthogonality of complex exponentials, we know that Then, we have where is the first kind Bessel function of order .
According to Debye’s relation [25], we know that decays exponentially when . From (13) we also know that when . Thus, we will have 2D Fourier representation of the HL condition as follows:
So far, we only assume that the object is supported on the unit disk. If the object is supported by a disk with a radius , then we replace by (where ) in (14) and (18) such that we can obtain the following conditions:where is the largest object support. Note that this property was also found in [18] by investigating the parallel-beam sinogram of a point object.
Figure 1 illustrates a double-wedge region of zero coefficients in the 2D Fourier transform of a sinogram, when .
(a)
(b)
(c)
2.4. HL Conditions in Fan-Beam Geometry
With variable substitution, that is, , ( is the rotation angle in fan-beam geometry, is the opening fan angle, and is the source-to-isocenter distance), the projection moment theorem (i.e., (5)) in fan-beam notation can be expressed aswhere
Similar to the parallel-beam case, we consider the object support is and expand in the Chebyshev-Fourier space as
Also, we readily obtain the relation , for Then, 2D Fourier transform of both sides of (23) and simplification according to (17) yield
According to (B.2) in Appendix B, the innermost integral with respect to becomes
It is known that the Weber-Schafheitlin’s integral decays very fast for [25]. Together with the relation , for , we finally arrive at
For an equally spaced fan-beam geometry, we apply the 2D Fourier transform to (23) with respect to and :where and denotes the source-detector distance.
Since is large compared to (for a small fan angle), we can make an approximation . Then, we will have a similar derivation as for (25):
Finally, we will obtain the HL conditions for an equally spaced fan-beam geometry:
Note that both (26) and (29) were approximately arrived at in [19] when investigating the fan-beam sinogram of a delta function point.
3. Data Extrapolation Using HL Consistency Conditions
The HL consistency conditions play an important role in image reconstruction as they can be used as a measure to restore consistency of the imperfect projection data (e.g., due to noise, motion, or truncation). In this work, we take advantage of the theory derived before for data extrapolation with truncated projection data. A flowchart of our proposed projection completion algorithm can be described as follows (also see Figure 2 for illustration).
Step 1. Segment the double-wedge region in the 2D Fourier transform of the sinogram as shown in Figure 1(c). Then, compute the slope of double wedge, that is, at double-wedge’s edge. The object support can be computed using (20); that is, , for instance, for a parallel-beam sinogram.
Step 2. Set up a shape model for sinogram completion (detruncation) and fit the model to the measured (truncated) data with a detruncation optimization algorithm that enforces the constraint that the values within the double-wedge region of the Fourier transformed sinogram are zero.(a)Set up a uniform ellipse model and initialize the model with two radii . For reasons of simplicity, assume a uniform density of the measured object. The density can be determined by two ways: heuristic preset, for example, water density; extrapolating the sinogram up to the object support, reconstructing it using an FBP framework, and using the mean of the reconstructed object as the density value. Additionally, current consistency conditions are only derived for a centered object; we assume the ellipse center .(b)Complete the truncated/measured sinogram using forward projections of the ellipse model. Here the forward projections only need to be computed for the regions outside of the scan’s field of view. where indicates the Radon Transform.(c)Optimize (adapt/smooth) the transition region of the completed sinogram by adding/subtracting an offset to each row of the forward-projected sinogram. The offset values are computed by comparing two neighboring pixels from the forward-projected and the measured sinogram as follows: where indicates the index of neighboring pixels between the extrapolated and the measured sinogram. Note that this rescaling/intensity adjustment of the forward-projected sinogram also weakens the impact of the ellipse density such that it is not mandatory to set the density as an additional optimization parameter.(d)Perform 2D fast Fourier transform (FFT) of the combined sinogram .(e)Sum up the energies in the double-wedge region where denotes the double-wedge region.(f)Minimize the energy in the double-wedge region by optimizing the ellipse parameters (g)Based on the estimated ellipse parameters , perform the sinogram completion and apply transition smoothing, as already described in (c).(h)Perform 2D FFT, set the residual energy in the double-wedge region to zero (hard constraint), and perform an inverse 2D FFT. To improve the reconstruction result, this step can be performed in an iterative manner.
Step 3. The resulting completed (detruncated) sinogram (being optimized in the Fourier domain) is reconstructed using any reconstruction algorithm that can be applied on nontruncated projection data, for example, the standard FBP framework.
The detection to the wedge region in Step 1 can be performed using edge-based segmentation with known properties such as line-symmetric for parallel-beam sinogram and point-symmetric for fan-beam case. For proof of concept, in current work we assume that object support is known and thus the wedge region can be directly computed using (29). Then, in (e) we set up a binary double-wedge mask with zero entries outside the double-wedge region and one within the double-wedge region. After applying the mask (element-wise multiplication) on the Fourier transformed completed sinogram from (d), we can sum up the energies in the double-wedge region.
For the task of minimization in (f), we use a Differential Evolution (DE) optimization [26] to search large spaces of candidate solutions and avoid the local minima. It is a stochastic, population-based global optimization method that appears fairly fast and robust for nondifferentiable and nonlinear objective functions. It uses a fixed number of parameter vectors as a population for each iteration (also referred to as a generation). Firstly, the trial parameter vectors (parent) are initialized on an interval which defines upper and lower bounds of parameters. At each iteration/generation, every parent is combined one by one with a set of new vectors which results in a set of trial vectors (children). These newly generated children vectors are mixed with a predetermined target vector with probability CR, generating the new trial vector. Finally, the new trial vector replaces the “least useful” parent that yields the largest cost function value if and only if its cost function value is lower than that of this parent. The steps above continue until some stopping criterion is reached. In this work, we choose Scheme DE 1 proposed in [26] with parametrization as , , and . We finish the optimization procedure if a preset maximum iteration number is reached.
4. Experiments and Results
The proposed method was validated and evaluated on the Shepp-Logan phantom and reprojected clinical data (data courtesy of CHI St. Luke’s Health-Baylor St. Luke’s Medical Center, Houston, TX, USA). All datasets were virtually collimated (by setting the area outside of the region of interest to zero) to a medium FOV and a small FOV. An equally spaced fan-beam imaging geometry is used. The sinogram of the Shepp-Logan phantom consists of 720 channels and 360 views over 360° rotation angle (full scan). Then, the phantom was reconstructed in a matrix with the radius of the medium ROI being 60 pixels and the radius of the small ROI being 30 pixels. The clinical data was forward-projected to a sinogram consisting of 1500 channels and 720 views over 360° rotation angle (full scan). The forward-projected sinogram was reconstructed in a matrix. The radius of the medium ROI is 95 pixels and the radius of the small ROI is 60 pixels. We also investigate the performance of the state-of-the-art water cylinder extrapolation method [2] and compare it with our proposed method.
Shown in Figures 3 and 4 are the reconstruction results from the Shepp-Logan phantom with different degrees of truncation (medium/small FOV). Difference images with respect to a reference FBP reconstruction from the original nontruncated data are also presented. The grayscale values are normalized into an arrangement of . As expected, reconstruction without any correction will generate typical cupping-like truncation artifacts and a substantial offset compared to the reference. These artifacts, on the other hand, are effectively compensated by using water cylinder extrapolation as well as the proposed method, yielding only a small error spreading over the difference images. Note that, for the Shepp-Logan phantom with a rather simple structure, the proposed method is able to accurately estimate the outline of the object, in contrast to inferior shape estimation from the water cylinder extrapolation scheme, in which extrapolation is performed by fulfilling the continuity assumptions. Figure 4 shows the case of severe truncation, where in general the reconstruction bias becomes larger as less data can be used for robust extrapolation and thus is more challenging for truncation correction algorithms. In this case, the reconstruction from the proposed method obviously shows less bias and truncation artifacts than that of the water cylinder extrapolation, while still retaining an accurate shape estimation. We also show the reconstruction images in a compressed window (see Figure 5) to demonstrate the robustness of the proposed method to both medium truncation and severe truncation case.
The results for reprojected clinical data are shown in Figures 6 and 7, also using two different degrees of truncation. The same reconstructions but with a compressed display window are shown in Figure 8. We can see that, for medium truncation (Figures 6(a) and 7(a)), both water cylinder extrapolation and the proposed method are able to reconstruct the image with high quality within the ROI: no typical cupping artifact and obvious bias are observed. However, when it comes to severe truncation, the water cylinder extrapolation performance degrades more than the performance of the proposed method. Figure 9 shows reconstructions of a reprojection from a slice of body scan, which basically have a similar observation to that of head scan. The line profiles in Figures 10 and 11 show that a bias always appears in water cylinder extrapolation results compared to the reference but is less observed in the proposed method. Quantitative results in Table 1 confirm that the performance of the water cylinder extrapolation appears to be ad hoc. It yields a root mean square error (RMSE) gray value of 271.6 while for the proposed method the error is always lower than 60 gray values. The proposed method is also able to nicely recover the structural information of the object and reduces the cupping artifacts. Thus, it yields a high correlation coefficient (CC) with respect to the reference, that is, 0.99, compared to a CC of 0.98 for the water cylinder extrapolation scheme.
(a)
(b)
(c)
5. Discussion
The Fourier-based Helgason-Ludwig consistency conditions are derived for both parallel-beam and fan-beam geometry, as described in (20), (26), and (29). The derivation outcomes indicate the Fourier property of a physically consistent sinogram: there are zero coefficients forming a double-wedge region (symmetric for parallel-beam and point-symmetric for fan-beam geometry) in its 2D Fourier transform. Interestingly, the same property was also observed previously in the literature by investigating the parallel-/fan-beam sinogram of a delta point object [18, 19]. In [19], the authors clarify that the approximation to a Bessel function (see Eq. in [19]) was arrived at intuitively and is validated empirically.
Motivated by these previous practical observations, we generate theoretical derivation. Our derivation stems directly from the original formulation of HL consistency conditions and is theoretically exact for parallel-beam geometry as well as equal-angle fan-beam geometry. In case of an equally spaced fan-beam case, such a derivation is not straightforward. To derive a similar property, we made an approximation that under the assumption that the opening fan angle is small. Therefore, there could potentially be a misestimation of zero-energy region for large fan angles, as also observed in [19].
The benefits of Fourier-based consistency conditions are as follows: First, rather than using a small subset of the consistency conditions as proposed in [2, 13–15], an infinite number of conditions are implicitly considered in the Fourier representation. Second, the 2D Fourier transform via FFT is computationally more efficient than other transforms (e.g., the Chebyshev-Fourier transform [16] or Lagrange-Fourier transform [11]). This allows us to develop more efficient sinogram recovery schemes, as demonstrated in this paper.
The sinogram-based extrapolation scheme we proposed in this work incorporates the Fourier-based consistency conditions as a constraint for optimizing the ellipse parameters so that the missing data can be more accurately fitted. Experiments on both phantom and clinical data yielded promising results. There are some limitations to this Fourier constrained extrapolation method applied to ROI reconstruction. First, the current derivation only involves the sinogram of a centered object. It is not clear how the zero energy region will change for off-center cases. Second, for evaluation we used a full scan fan-beam geometry. We observed that, for a short scan acquisition where projection data is acquired only over a range of plus the fan angle, some nonzero values also appear in zero energy region, which may affect the optimization procedure. Thus, corresponding consistency conditions that also account for an off-centered object and short scan acquisition would be interesting for future work. Also, in this work we used a uniform ellipse model to generate the sinogram outside measured region. Such an assumption is well suited for the imaged objects that can be approximated by a single ellipse, for example, head scan. For complicated objects such as knee scan, multiple ellipses may be superposed and optimized.
6. Conclusion
In this paper, we theoretically derived the 2D Fourier-based Helgason-Ludwig consistency conditions that can be evaluated very efficiently via FFT. Then, we suggested a sinogram-based extrapolation scheme that incorporates these consistency conditions as a constraint for optimizing the ellipse parameters. Then, the forward projection of the optimized ellipse can be used to complete the truncation data. Experiments on both phantom and clinical data yielded promising results of the proposed approach. The reconstruction results indicate that the proposed approach substantially outperforms a conventional water cylinder extrapolation approach [2], particularly for severe truncation, regarding both image quality and residual artifacts.
Appendix
A. Proving That Form an Orthogonal Basis of
B. Relationship between Bessel Function and Chebyshev Polynomial
The Fourier transform of the Bessel function can be computed as for . We can swap the variables and and move the term to the left side: where
Then, substituting by and computing the Fourier transform of both sides give
Disclosure
The concepts and information presented in this paper are based on research and are not commercially available.
Conflicts of Interest
Martin Berger, Shiyang Hu, Andre Aichert, and Andreas Maier have no conflicts of interest. Yan Xia was supported by Siemens Healthcare GmbH, Forchheim, Germany, and by Erlangen Graduate School in Advanced Optical Technologies (SAOT). Sebastian Bauer is an employee of Siemens Healthcare GmbH, Forchheim, Germany.
Acknowledgments
The authors gratefully acknowledge funding by Siemens Healthcare GmbH and of the Erlangen Graduate School in Advanced Optical Technologies (SAOT) by the German Research Foundation (DFG) in the framework of the German Excellence Initiative.