Derivative-Hilbert-Backprojection based image reconstruction from truncated projections in helical cone-beam CT

NUCLEAR ELECTRONICS AND INSTRUMENTATION

Derivative-Hilbert-Backprojection based image reconstruction from truncated projections in helical cone-beam CT

ZHANG Feng，

YAN Bin，

LI Lei，

XI Xiao-Qi，

JIANG Hua，

WEI Xing，

ZHANG Xiang，

CUI Jin-Xian

Nuclear Science and Techniques

Vol.26, No.2

Article number 020401

Published in print 20 Apr 2015

Available online 20 Apr 2015

DOI：10.13538/j.1001-8042/nst.26.020401

42804

In helical cone-beam computed tomography (CT), Feldkamp-Davis-Kress (FDK) based image reconstruction algorithms are by far the most popular. However, artifacts are commonly met in the presence of lateral projection truncation. The reason is that the ramp filter is global. To restrain the truncation artifacts, an approximate reconstruction formula is proposed based on the Derivative-Hilbert-Backprojection (DHB) framework. In the method, the first order derivative filter is followed by the Hilbert transform. Since the filtered projection values are almost zero by the first order derivative filter, the following Hilbert transform has little influence on the projection values, even though the projections are laterally truncated. The proposed method has two main advantages. First, it has comparable computational efficiency and image quality as well as the conventional helical FDK algorithm for non-truncated projections. The second advantage is that images can be reconstructed with acceptable quality and much lower computational cost in comparison to the Laplace operator based algorithm in cases with truncated projections. To point out the advantages of our method, simulations on the computer and real data experiments on our laboratory industrial cone-beam CT are conducted. The simulated and experimental results demonstrate that the method is feasible for image reconstruction in the case of projection truncation.

Helical cone-beam CTImage reconstructionTruncation artifactsDerivative-Hilbert-Backprojection

I. INTRODUCTION

In the field of cone-beam computed tomography (CT), helical scans have the virtue of higher scan speed and more uniform axial errors compared to the circular scan [1]. The image reconstruction problem of the helical scan has been well solved by means of approximate or exact algorithms [2-4]. Among these image reconstruction methods, Wang et al.’s method [2] (called the helical FDK method hereafter) is computationally efficient and widely used in practice. However, the method is only suitable in cases that the objects under investigation are surrounded by the field of view (FOV) defined by the detector size, scan magnification and rotation axis location. In other words, when the object exceeds the FOV, the sampled projections are truncated and the method becomes challenged. When the projections are truncated on both sides, called the interior problem, the image reconstruction has been proven to have no unique or unstable solutions [5]. If the helical FDK method was employed to deal with the truncated projections, the resulting reconstruction image quality will be severely degraded by truncation artifacts [6].

With the development of a local reconstruction algorithm, there are new solutions to the interior problem [7, 8, 9, 10, 11, 12, 13, 14, 15]. Of these solutions, the approximate truncation reconstruction computed tomography (ATRACT)-type algorithms researched by Dennerlein et al.’s group are highly practical [10, 11, 12, 13, 14, 15]. Its main idea is to use a cascade of local and global filters to replace the ramp filter in the FDK-type algorithms [16, 17]. Thus, this kind of algorithm is expected to have the same advantages of high computational efficiency and low resource demand as those of the FDK-type methods. Based on Dennerlein et al.’s work, Zou et al. [4] proposed a general Laplace operator based reconstruction algorithm (LORA) for helical scans recently (called the Laplace operator based method in this paper), where the Laplace operator and 2D Radon based filters are regarded as the local and global filters, respectively. Although the method can provide better reconstruction over the helical FDK algorithm for a truncated spiral cone beam CT, the Radon based filter is quite time-consuming. In fact, many improvements have been made on the computational efficiency since the publication of Dennerlein et al.’s first work [11, 13, 15, 19]. Among these improvements researches, Wang et al. [19] proposed to filter the projection by the first order derivative and Hilbert filters, which can be described as the Derivative-Hilbert-Backprojection (DHB) framework. In terms of computational cost and image quality, Wang et al.’s method is outstanding.

In this paper, a novel approximate reconstruction formula based on Wang et al.’s result is proposed for truncated projections of helical cone-beam CT. The remainder of the paper is organized as follows. The geometry of a helical scan is described in section II. The theory of our method is introduced in detail in section III. Section IV discusses the numerical simulation and experimental results. Section V is Conclusion.

II. NOTATION AND BASIC DEFINITIONS

The helical scan geometry is illustrated in Fig. 1. Without loss of generality, the X-ray source and the flat panel detector are assumed to rotate around the imaging object in the data acquisition of a helical scan. The actual position of the source, which is parameterized by the azimuth angle $λ \in ℝ$ , where $ℝ$ is the real number space, is denoted by ${\vec{r}}_{0} (λ)$ . Three coordinate systems are hereby included in the geometry. The first one is the 3-dimensional (3D) Cartesian coordinate system, Oxyz, which is fixed in the helical scan. The second one is the rotational coordinate system which is spanned by ${\vec{e}}_{u} (λ) = (-,,0)$ , ${\vec{e}}_{ν} (λ = (0,0,1)$ and ${\vec{e}}_{w} (λ) = (,,0)$ . The third one is the 2-dimensional (2D) Cartesian coordinate system on the flat detector with u- and ν- axes. For a helical trajectory with a constant pitch, h, and constant radius, R, one has ${\vec{r}}_{0} (λ) = (R, R, h λ / (2 π))$ . The pitch is defined as the distance in the z-direction of two points on the helix that are exactly one turn apart. The flat panel detector is arranged to be parallel to the vectors ${\vec{e}}_{u} (λ)$ and ${\vec{e}}_{ν} (λ)$ , and orthogonal to the vector ${\vec{e}}_{w} (λ)$ . The distance from the source to the flat panel detector is D. R represents the distance between the source and the rotation axis. The FOV of helical cone-beam CT is a bounding cylinder C.

Fig. 1.

The geometry of helical cone-beam CT.

It is assumed that the spatial distribution of the X-ray absorption coefficient of imaging object is a compact function denoted by $f (\vec{r})$ , which is bounded by Ω. The projection of object $f (\vec{r})$ onto the flat panel detector is

g (λ, u, ν) = \int_{R} (f ({\vec{r}}_{0} (λ) + \vec{θ} t)) d t .

(1)

The direction of the X-ray given by the parameters λ, u and ν is defined by

\vec{θ} = \frac{u {\vec{e}}_{u} (λ) + ν {\vec{e}}_{ν} (λ) - D {\vec{e}}_{w} (λ)}{\sqrt{u^{2} + ν^{2} + D^{2}}},

(2)

where the coordinate (u,ν) on the intersection where the X-ray passing through $\vec{r} = (x, y, z)$ hits the detector can be computed with

u = \frac{\vec{r} \cdot {\vec{e}}_{u} (λ)}{R - \vec{r} \cdot {\vec{e}}_{w} (λ)} D and ν = \frac{\vec{r} \cdot {\vec{e}}_{ν} (λ)}{R - \vec{r} \cdot {\vec{e}}_{w} (λ)} D .

(3)

If the object $f (\vec{r})$ cannot be laterally surrounded by the FOV of a helical scan, i.e. Ω cannot be contained by C, the captured projection is laterally truncated.

III. THEORY

To obtain better image quality from transversely truncated projections in a helical trajectory, a novel reconstruction algorithm is proposed based on Wang et al.’s work. The new reconstruction formula can be written as

f (\vec{r}) = \frac{1}{4 π} \int_{θ - π}^{θ + π} \frac{R D}{{[R - \vec{r} \cdot {\vec{e}}_{w} (λ)]}^{2}} [\int_{- \infty}^{\infty} h_{H} (u - u^{'}) \frac{\partial}{\partial u} [\frac{D}{\sqrt{D^{2} + u^{2} + ν^{2}}} g (λ, u, ν)] d u^{'}] d λ,

(4)

where θ = 2πz/h. Formula (4) can be implemented by the following four steps.

Step 1. Weighing the projection data

g_{1} (λ, u, ν) = \frac{D}{\sqrt{D^{2} + u^{2} + ν^{2}}} g (λ, u, ν),

(5)

where g(λ, u, ν) is the projection at a view of λ and (u,ν) is a coordinate of a detector pixel.

Step 2. First order derivation of the weighted projection with respect to u

g_{2} (λ, u, ν) = \frac{\partial}{\partial u} g_{1} (λ, u, ν) .

(6)

Step 3. Hilbert filtering of g₂(λ, u, ν) row by row

g_{3} (λ, u, ν) = \frac{1}{2 π} g_{2} (λ, u, ν) * h_{H} (u),

(7)

where hH = 1/(πu) is the Hilbert transform kernel and the symbol * signifies a convolution operation.

Step 4. Back-projecting the filtered projections onto the 3D volume

f (\vec{r}) = \frac{1}{2} \int_{θ - π}^{θ + π} \frac{R D}{{[R - \vec{r} \cdot {\vec{e}}_{w} (λ)]}^{2}} g_{3} (λ, u, ν) d λ .

(8)

Mathematically, formula (4) is equivalent to formula (10) in reference [2]. In the helical FDK algorithm, the filtration kernel |ω| is ramp-like in the frequency domain.

In our method, the ramp filter is replaced by a cascade of first order derivative and Hilbert operators. Hence, the proposed reconstruction method has the same characteristics as the helical FDK algorithm for global reconstruction. To obtain the correct filtered result, the acquired projection must be complete because |ω| is global. When the projections are transversely truncated, the filtered values are not equal to those from non-truncated projections. As a result, truncation artifacts resulting from projection truncation arise in the reconstruction image. However, the filtered projection values, g₂(λ,u,ν), by the first order derivative operator are close to zeroes in most regions. Even if the following Hilbert filter is still global, the final filtered errors can be significantly reduced in comparison to the helical FDK algorithm.

Different from our method, Zou et al.’ method utilized the Laplace operator and 2D Radon based filters [18]. When the weighted projection $g_{1}^{Z o u} (λ, u, ν)$ is available, Zou et al. performed the Laplace operator based filtering by

g_{2}^{Z o u} (λ, u, ν) = [\frac{\partial^{2}}{\partial u^{2}} + \frac{\partial^{2}}{\partial ν^{2}}] g_{1}^{Z o u} (λ, u, ν),

(9)

and 2D Radon based filtering by

g_{F}^{Z o u} (λ, u, ν) = - \frac{1}{4 π^{2}} \int_{0}^{π} | \cos μ | g_{3}^{Z o u} (λ, u, ν) d μ,

(10)

where

g_{3}^{Z o u} (λ, u^{*}, ν^{*}) = \int_{Ω_{λ}} g_{2}^{Z o u} (λ, u, ν) δ (u \cos μ + ν \sin μ - s) d u d μ,

(11)

is the Radon transform, in which s = u^* cos μ + ν^* sin μ. In the implementation of the 2D Radon based filter, discretions of s and λ (denoted by Δs and Δλ, respectively) are critical to the quality of the filtered projection. Obviously, the smaller the values of Δs and Δλ are, the better the quality of filtered projection can be achieved. Meanwhile, more computation time is also demanded. In our following discussions, Δs=0.5 pixel and Δλ=0.5°.

To compare the filtered results by the three reconstruction methods mentioned above, the non-truncated 2D Shepp-Logan phantom with 512x512 pixels and the truncated one with 256x512 pixels are generated and filtered, respectively. The filtered results of the non-truncated phantom image by the three methods are shown in Fig. 2. From the results, we see that the filtered image using Zou et al.’s method has a jagged shape. In comparison, the filtered result using our method basically coincides with that using the ramp filter in the helical FDK algorithm.

Fig. 2.

The upper-left image (a) is the non-truncated 2D Shepp-Logan phantom. The following three images [(b),(c),(d)] are filtered results using the helical FDK, Zou et al.’s method and the proposed method. The display window of the filtered results is [–0.04, 0.05]. The rightmost image (e) illustrates the profiles of the middle row data of the three filtered images.

Also, the truncated 2D Shepp-Logan phantom is filtered by the three kinds of filters and the results are shown in Fig. 3. Because the phantom image is not complete, the filtered image using the helical FDK algorithm deviates from the true value especially at the truncated places. Like the results in the first simulation, the filtered result using Zou et al.’s method has a jagged shape. Meanwhile, the result using our method coincides well with the true value.

Fig. 3.

The leftmost image is the truncated 2D Shepp-Logan phantom. The following three images are filtered results using the helical FDK, Zou et al.’s method and the proposed method. The display window of the filtered results is [–0.04, 0.05]. The rightmost image illustrates the profiles of the middle row data of the three filtered images.

It is worth noting that filtering one image with 512x512 pixels by Zou et al.’s method consumes 25.609 s and 15.703 s for one image with 256x512 pixels. In contrast, our method consumes 0.204 seconds and 0.109 seconds respectively, and the helical FDK consumes 0.094 s and 0.046 s. From the results above, we believe that our method has the advantages of less time cost over Zou et al.’s method and higher quality over the helical FDK algorithm for truncated projections.

In fact, the DHB based reconstruction method was proposed long ago. At that time, the first order derivative of DHB was thought to be sensitive to noise and harmful to the reconstruction image quality. Therefore, researchers and engineers preferred the Filtration Back-projection (FBP) to the DHB method. To implement the DHB based method in practice, projection preprocessing such as projection denoising is conducted to restrain the noise of projection before image reconstruction in our implementation.

IV. EXPERIMENTS AND DISCUSSIONS

To verify the efficiency of the proposed method, numerical simulations and real data experiments were conducted, where the helical FDK algorithm, Zou et al.’s method, and our method were all employed to deal with the acquired projections. These algorithms were implemented on a HP workstation with 32GB memory, Intel E5-2620 CPU, and programmed with C++ language on Microsoft Visual Studio 2010.

A. Simulated data

In the following experiment, we demonstrated the efficiency of our method on simulated projection data from the 3D Shepp-Logan phantom [20]. The 3D Shepp-Logan phantom size was 25.6 mm×25.6 mm×25.6 mm. Image reconstructions from non-truncated (Simulation 1) and truncated projections (Simulation 2) are discussed, respectively. In the Simulation 1, non-truncated projections along a helical trajectory were generated, in which the detailed simulation parameters are shown in the second column of TABLE 1. In Simulation 2, truncated projections along the same helical trajectory were generated and the simulation parameters are illustrated in the third column of TABLE 1. Because the width of the detector is 256 pixels in the Simulation 2, the simulated projections are laterally truncated.

Parameters of numerical simulations

parameters	Simulation 1	Simulation 2
source to	150 mm	150 mm
object distance
source to	600 mm	600 mm
detector distance
size of detector	256 pixel×256 pixel	256 pixel×512 pixel
size of detector cell	0.2 mm×0.2 mm	0.2 mm×0.2 mm
number of	360	360
projections per turn
pitch of helix	25 mm	25 mm
Number of turn	2	2

In the Simulation 1, the reconstruction matrix was 512x512x512 and the voxel edge length was 0.05 mm along each axis. The reconstruction images from the simulated non-truncated projections are shown in Fig. 4. Fig. 5 illustrates the profiles of middle row data of the first row in Fig. 4. From the results, it was found that the three methods can provide high quality reconstruction images for non-truncated projections. Because the filtered projection is rough by the Laplace operator based method, the profile of the middle row data is jagged. By comparison, reconstruction image quality using our method is slightly superior to that of Zou et al.’s method. To analyze the reconstruction error quantitatively, the root mean square errors (RMSEs) of the three results relative to the phantom were calculated. The resulting RMSEs were 1.128×10^-4, 1.149×10^-4 and 1.147×10^-4, respectively. Additionally, the consumed times for reconstructing a 512x512 pixel slice by the three methods was recorded and is shown in TABLE 2. We can see that Zou et al.’s method is the most time consuming and that the time for our method is slightly longer than the helical FDK algorithm. The reason is that 2D Radon based filtering is very computationally demanding, which has been stated in Section III.

Fig. 4.

The images in the first column are slices of the phantom. The following column images from left to right are reconstructed from non-truncated projections by use of the helical FDK algorithm, Zou et al.’s method and the proposed method, respectively. The images from up to down are slices at z=12.8 mm, y=12.8 mm and x=12.8 mm, respectively. The display window is [0, 1].

Fig. 5.

The top image illustrates the profiles of middle row data of slices at z=12.8 mm. The bottom image illustrates the profiles of middle column data of slices at x=12.8 mm.

Calculation times of the three reconstruction algorithms in seconds

Simulation	The helical FDK algorithm	LORA method	The proposed method
Simulation 1	64	9232	92
Simulation 2	38	5666	52

In Simulation 2, the reconstruction matrix was 256x256x512 and the voxel edge length was 0.05 mm along each axis. During the filtering process the mean value of the truncated projection data was removed by Zou et al.’s method and our method, as a result, an offset artifact or bias-like artifact arose. For convenient comparison, the reconstruction images using the two methods were corrected by scaling and offset corrections. The two correction techniques can be found in reference [15]. The reconstruction images from the simulated truncated projections are shown in Fig. 6. From the results, one can find that the reconstruction image using the helical FDK algorithm is severely contaminated by truncation artifacts. Since the filtered projection values using the Laplace operator based filter are closer to zero than those using the first order derivative filter, the reconstruction image quality on the edge is slightly better. To evaluate on the quality of the reconstruction images quantitatively, signal-to-noise ratio (SNR) was computed. The region of interest (ROI) is depicted by the white circle in Fig. 6. The SNRs of ROIs by the helical FDK algorithm, Zou et al.’s method, and our method were 10.38, 19.54 and 21.18 dB, respectively. We also recorded the computation time for reconstructing a 256x256 pixel slice and the results are shown in TABLE 2. The conclusions are similar to Simulation 1.

Fig. 6.

The left image in the first row is the slice of phantom. The following images are reconstructed from truncated projections by use of the helical FDK algorithm, Zou et al.’s method and the proposed method, respectively. The display window is [0, 1]. The last image illustrates the profiles of the middle row data of the three reconstruction images.

B. Real data

Real data experiments were conducted on our laboratory industrial cone-beam CT, which consists of a micro-focus X-ray source (FXT-225.48, Yxlon, Germany), a high-precision manipulator and a flat panel detector with 2304x3200 pixels (4030E, Varian, America). A car spark plug was scanned along a helical trajectory. In the experiment, the source to object distance was set to 125 mm and the source to detector distance was 1000 mm, respectively. To obtain high quality reconstruction, the pitch of the helical scan was set to 36 mm. 360 projections with a resolution of 0.127 mm/pixel from the car spark was sampled in one turn. The turn number of helical scan is 3, so 1080 projections are collected. To avoid the geometry artifacts interference, the sampled projections are firstly corrected using a practical geometric calibration method [1]. It is worth noting that the collected projections are non-truncated since the true values from car spark plug are unknown. For convenient comparison, the helical FDK algorithm was employed to the reconstruct the image from non-truncated projections and the reconstruction image was regarded as the true value. The reconstructed 3D volume consists of 1536x1536x5120 voxels with 0.018×0.018×0.018 mm per voxel. A photo of the car spark plug and the reconstruction images are shown in Fig. 7.

Fig. 7.

The left upper image is the photo of the car spark plug. The right upper one is the 768 slice of the 3D volume along y-axis. The bottom one is the 2400 slice along z-axis depicted by the white line in the middle image. The display window is [0, 1].

The performance of the three methods was evaluated on truncated projections. To obtain the truncated projections, the central 512 columns of all projections were cropped virtually. Just as in the simulation process, the helical FDK algorithm, Zou et al.’s method, and the proposed method were all employed to deal with the cropped projections. The reconstruction matrix was 512x512x5120 and the voxel size was 0.018×0.018×0.018 mm. The slices at y=256 and z=2400 are shown in Figs. 8 and 9, respectively. From the results, we see that the reconstruction image from truncated projections using the helical FDK algorithm was severely influenced by truncation artifacts. However, Zou et al.’s method and our method can significantly reduce the truncation artifacts and provide comparable image quality. Also, the SNRs of ROIs using these three methods were 9.81, 14.06 and 23.06 dB, respectively. The ROI is depicted by the white circle in Fig. 9. Similar to the conclusions in the simulations, our method was still faster than Zou et al.’s method. For the sake of brevity, comparison of the computational time is omitted here.

Fig. 8.

The top image is cut out from the middle image in Fig. 4. The following images from up to down are reconstructed from truncated projections by means of the helical FDK algorithm, Zou et al.’s method and the proposed method.

Fig. 9.

The left image in the upper row is cut out from the rightmost image in Fig. 4. The following images are reconstructed by use of the helical FDK algorithm, Zou et al.’s method and the proposed method, respectively. The display window is [0, 1]. The last image in the lower row illustrates the profiles of the middle row data of the images above.

V. CONCLUSION

For global tomography, the helical FDK algorithm is by far the most popular one in helical scan reconstructions. However, the reconstruction image quality is always degraded by truncation artifacts in the presence of projection truncation. In this paper, a novel approximate reconstruction formula based upon DHB framework was proposed for helical interior cone-beam CT. The method performed the ramp filtering operation in the helical FDK algorithm by applying a cascade of first order derivative filtering and Hilbert filtering. Because the intermediate filtered projection values using the first order derivative filter were close to zero, the final filtered results were not significantly changed by the following Hilbert filter. Theoretically, the method is identical to the helical FDK algorithm for non-truncated cases. Accordingly, our method has the same advantages, such as low resource demand, as the helical FDK algorithm.

We also compared the method with the Laplace operator based method proposed by Zou et al.. The experimental results showed that our method possessed the advantages of less time cost and comparable image quality. Hence, our method has a significant practical implication. Since the first order derivative operation was included in our method, it may be sensitive to noise in practice. However, this potential weakness can be overcome by introducing a projection denoising in practice, such as non-local means.

References

[1]

F Zhang, B Yan, L Li, et al.

Practical geometric calibration for helical cone-beam industrial computed tomography

. J X-ray Sci Technol, 2014, 22: 19-35. DOI: 10.3233/XST-130406