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

Derivative-Hilbert-Backprojection based image reconstruction from truncated projections in helical cone-beam CT 1 first-author 1 corresp Corresponding author, tom.yan@gmail.com Corresponding author, tom.yan@gmail.com tom.yan@gmail.com 1 1 1 1 1 1 National Digital Switching System Engineering & Technology Research Center, Zhengzhou 450002, China National Digital Switching System Engineering & Technology Research Center, Zhengzhou 450002, China 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 CT Image reconstruction Truncation artifacts Derivative-Hilbert-Backprojection I. INTRODUCTION 1 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 1 ]. The image reconstruction problem of the helical scan has been well solved by means of approximate or exact algorithms [ 2 2 - 4 4 2-4 ]. Among these image reconstruction methods, Wang et al .’s method [ 2 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 vi ew (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 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 6 ]. With the development of a local reconstruction algorithm, there are new solutions to the interior problem [ 7 7 , 8 8 , 9 9 , 10 10 , 11 11 , 12 12 , 13 13 , 14 14 , 15 15 ]. Of these solutions, the approximate truncation reconstruction computed tomography (ATRACT)-type algorithms researched by Dennerlein et al .’s group are highly practical [ 10 10 , 11 11 , 12 12 , 13 13 , 14 14 , 15 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 16 , 17 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 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 11 , 13 13 , 15 15 , 19 19 ]. Among these improvements researches, Wang et al . [ 19 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 1 The helical scan geometry is illustrated in Fig. 1 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 <math id="M1"><mrow><mi>λ</mi><mo>∈</mo><mi>ℝ</mi></mrow></math> , where <math id="M2"><mi>ℝ</mi></math> is the real number space, is denoted by <math id="M3"><mrow><msub><mover accent="true"><mi>r</mi><mo>→</mo></mover><mn>0</mn></msub><mo stretchy="false">(</mo><mi>λ</mi><mo stretchy="false">)</mo></mrow></math> . Three coordinate systems are hereby included in the geometry. The first one is the 3-dimensional (3D) Cartesian coordinate system, O xyz , which is fixed in the helical scan. The second one is the rotational coordinate system which is spanned by <math id="M4"><mrow><msub><mover accent="true"><mi>e</mi><mo>→</mo></mover><mi>u</mi></msub><mo stretchy="false">(</mo><mi>λ</mi><mo stretchy="false">)</mo><mo>=</mo><mo stretchy="false">(</mo><mo>−</mo><mn>,,0</mn><mo stretchy="false">)</mo></mrow></math> , <math id="M5"><mrow><msub><mover accent="true"><mi>e</mi><mo>→</mo></mover><mi>ν</mi></msub><mo stretchy="false">(</mo><mi>λ</mi><mo>=</mo><mo stretchy="false">(</mo><mn>0,0,1</mn><mo stretchy="false">)</mo></mrow></math> and <math id="M6"><mrow><msub><mover accent="true"><mi>e</mi><mo>→</mo></mover><mi>w</mi></msub><mo stretchy="false">(</mo><mi>λ</mi><mo stretchy="false">)</mo><mo>=</mo><mo stretchy="false">(</mo><mn>,,0</mn><mo stretchy="false">)</mo></mrow></math> . 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 <math id="M7"><mrow><msub><mover accent="true"><mi>r</mi><mo>→</mo></mover><mn>0</mn></msub><mo stretchy="false">(</mo><mi>λ</mi><mo stretchy="false">)</mo><mo>=</mo><mo stretchy="false">(</mo><mi>R</mi><mn>,</mn><mi>R</mi><mn>,</mn><mi>h</mi><mi>λ</mi><mo>/</mo><mo stretchy="false">(</mo><mn>2</mn><mi>π</mi><mo stretchy="false">)</mo><mo stretchy="false">)</mo></mrow></math> . 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 <math id="M8"><mrow><msub><mover accent="true"><mi>e</mi><mo>→</mo></mover><mi>u</mi></msub><mo stretchy="false">(</mo><mi>λ</mi><mo stretchy="false">)</mo></mrow></math> and <math id="M9"><mrow><msub><mover accent="true"><mi>e</mi><mo>→</mo></mover><mi>ν</mi></msub><mo stretchy="false">(</mo><mi>λ</mi><mo stretchy="false">)</mo></mrow></math> , and orthogonal to the vector <math id="M10"><mrow><msub><mover accent="true"><mi>e</mi><mo>→</mo></mover><mi>w</mi></msub><mo stretchy="false">(</mo><mi>λ</mi><mo stretchy="false">)</mo></mrow></math> . 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 . It is assumed that the spatial distribution of the X-ray absorption coefficient of imaging object is a compact function denoted by <math id="M11"><mrow><mi>f</mi><mo stretchy="false">(</mo><mover accent="true"><mi>r</mi><mo>→</mo></mover><mo stretchy="false">)</mo></mrow></math> , which is bounded by Ω . The projection of object <math id="M12"><mrow><mi>f</mi><mo stretchy="false">(</mo><mover accent="true"><mi>r</mi><mo>→</mo></mover><mo stretchy="false">)</mo></mrow></math> onto the flat panel detector is <math display="block" id="M13"><mrow><mi>g</mi><mo stretchy="false">(</mo><mi>λ</mi><mn>,</mn><mi>u</mi><mn>,</mn><mi>ν</mi><mo stretchy="false">)</mo><mo>=</mo><mstyle displaystyle="true"><mrow><msub><mo>∫</mo><mi>R</mi></msub></mrow></mstyle><mrow><mo>(</mo><mrow><mi>f</mi><mo stretchy="false">(</mo><msub><mover accent="true"><mi>r</mi><mo>→</mo></mover><mn>0</mn></msub><mo stretchy="false">(</mo><mi>λ</mi><mo stretchy="false">)</mo><mo>+</mo><mover accent="true"><mi>θ</mi><mo>→</mo></mover><mi>t</mi><mo stretchy="false">)</mo></mrow><mo>)</mo></mrow><mtext>d</mtext><mi>t</mi><mn>.</mn></mrow></math> The direction of the X-ray given by the parameters λ , u and ν is defined by <math display="block" id="M14"><mrow><mover accent="true"><mi>θ</mi><mo>→</mo></mover><mo>=</mo><mfrac><mrow><mi>u</mi><msub><mover accent="true"><mi>e</mi><mo>→</mo></mover><mi>u</mi></msub><mo stretchy="false">(</mo><mi>λ</mi><mo stretchy="false">)</mo><mo>+</mo><mi>ν</mi><msub><mover accent="true"><mi>e</mi><mo>→</mo></mover><mi>ν</mi></msub><mo stretchy="false">(</mo><mi>λ</mi><mo stretchy="false">)</mo><mo>−</mo><mi>D</mi><msub><mover accent="true"><mi>e</mi><mo>→</mo></mover><mi>w</mi></msub><mo stretchy="false">(</mo><mi>λ</mi><mo stretchy="false">)</mo></mrow><mrow><msqrt><mrow><msup><mi>u</mi><mn>2</mn></msup><mo>+</mo><msup><mi>ν</mi><mn>2</mn></msup><mo>+</mo><msup><mi>D</mi><mn>2</mn></msup></mrow></msqrt></mrow></mfrac><mn>,</mn></mrow></math> where the coordinate ( u , ν ) on the intersection where the X-ray passing through <math id="M15"><mrow><mover accent="true"><mi>r</mi><mo>→</mo></mover><mo>=</mo><mo stretchy="false">(</mo><mi>x</mi><mn>,</mn><mi>y</mi><mn>,</mn><mi>z</mi><mo stretchy="false">)</mo></mrow></math> hits the detector can be computed with <math display="block" id="M16"><mrow><mi>u</mi><mo>=</mo><mfrac><mrow><mover accent="true"><mi>r</mi><mo>→</mo></mover><mo>⋅</mo><msub><mover accent="true"><mi>e</mi><mo>→</mo></mover><mi>u</mi></msub><mo stretchy="false">(</mo><mi>λ</mi><mo stretchy="false">)</mo></mrow><mrow><mi>R</mi><mo>−</mo><mover accent="true"><mi>r</mi><mo>→</mo></mover><mo>⋅</mo><msub><mover accent="true"><mi>e</mi><mo>→</mo></mover><mi>w</mi></msub><mo stretchy="false">(</mo><mi>λ</mi><mo stretchy="false">)</mo></mrow></mfrac><mi>D</mi><mtext> </mtext><mtext>and</mtext><mtext> </mtext><mi>ν</mi><mo>=</mo><mfrac><mrow><mover accent="true"><mi>r</mi><mo>→</mo></mover><mo>⋅</mo><msub><mover accent="true"><mi>e</mi><mo>→</mo></mover><mi>ν</mi></msub><mo stretchy="false">(</mo><mi>λ</mi><mo stretchy="false">)</mo></mrow><mrow><mi>R</mi><mo>−</mo><mover accent="true"><mi>r</mi><mo>→</mo></mover><mo>⋅</mo><msub><mover accent="true"><mi>e</mi><mo>→</mo></mover><mi>w</mi></msub><mo stretchy="false">(</mo><mi>λ</mi><mo stretchy="false">)</mo></mrow></mfrac><mi>D</mi><mn>.</mn></mrow></math> If the object <math id="M17"><mrow><mi>f</mi><mo stretchy="false">(</mo><mover accent="true"><mi>r</mi><mo>→</mo></mover><mo stretchy="false">)</mo></mrow></math> 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 1 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 <math display="block" id="M18"><mrow><mi>f</mi><mo stretchy="false">(</mo><mover accent="true"><mi>r</mi><mo>→</mo></mover><mo stretchy="false">)</mo><mo>=</mo><mfrac><mn>1</mn><mrow><mn>4</mn><mi>π</mi></mrow></mfrac><mstyle displaystyle="true"><mrow><msubsup><mo>∫</mo><mrow><mi>θ</mi><mo>−</mo><mi>π</mi></mrow><mrow><mi>θ</mi><mo>+</mo><mi>π</mi></mrow></msubsup></mrow></mstyle><mfrac><mrow><mi>R</mi><mi>D</mi></mrow><mrow><msup><mrow><mrow><mo>[</mo><mrow><mi>R</mi><mo>−</mo><mover accent="true"><mi>r</mi><mo>→</mo></mover><mo>⋅</mo><msub><mover accent="true"><mi>e</mi><mo>→</mo></mover><mi>w</mi></msub><mo stretchy="false">(</mo><mi>λ</mi><mo stretchy="false">)</mo></mrow><mo>]</mo></mrow></mrow><mn>2</mn></msup></mrow></mfrac><mrow><mo>[</mo><mrow><mstyle displaystyle="true"><mrow><msubsup><mo>∫</mo><mrow><mo>−</mo><mi>∞</mi></mrow><mi>∞</mi></msubsup></mrow></mstyle><msub><mi>h</mi><mi>H</mi></msub><mo stretchy="false">(</mo><mi>u</mi><mo>−</mo><msup><mi>u</mi><mo>′</mo></msup><mo stretchy="false">)</mo><mfrac><mo>∂</mo><mrow><mo>∂</mo><mi>u</mi></mrow></mfrac><mrow><mo>[</mo><mrow><mfrac><mi>D</mi><mrow><msqrt><mrow><msup><mi>D</mi><mn>2</mn></msup><mo>+</mo><msup><mi>u</mi><mn>2</mn></msup><mo>+</mo><msup><mi>ν</mi><mn>2</mn></msup></mrow></msqrt></mrow></mfrac><mi>g</mi><mo stretchy="false">(</mo><mi>λ</mi><mn>,</mn><mi>u</mi><mn>,</mn><mi>ν</mi><mo stretchy="false">)</mo></mrow><mo>]</mo></mrow><mtext>d</mtext><msup><mi>u</mi><mo>′</mo></msup></mrow><mo>]</mo></mrow><mtext>d</mtext><mi>λ</mi><mn>,</mn></mrow></math> where θ = 2 πz/h . Formula (4) can be implemented by the following four steps. Step 1. Weighing the projection data <math display="block" id="M19"><mrow><msub><mi>g</mi><mn>1</mn></msub><mo stretchy="false">(</mo><mi>λ</mi><mn>,</mn><mi>u</mi><mn>,</mn><mi>ν</mi><mo stretchy="false">)</mo><mo>=</mo><mfrac><mi>D</mi><mrow><msqrt><mrow><msup><mi>D</mi><mn>2</mn></msup><mo>+</mo><msup><mi>u</mi><mn>2</mn></msup><mo>+</mo><msup><mi>ν</mi><mn>2</mn></msup></mrow></msqrt></mrow></mfrac><mi>g</mi><mo stretchy="false">(</mo><mi>λ</mi><mn>,</mn><mi>u</mi><mn>,</mn><mi>ν</mi><mo stretchy="false">)</mo><mn>,</mn></mrow></math> 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 <math display="block" id="M20"><mrow><msub><mi>g</mi><mn>2</mn></msub><mo stretchy="false">(</mo><mi>λ</mi><mn>,</mn><mi>u</mi><mn>,</mn><mi>ν</mi><mo stretchy="false">)</mo><mo>=</mo><mfrac><mo>∂</mo><mrow><mo>∂</mo><mi>u</mi></mrow></mfrac><msub><mi>g</mi><mn>1</mn></msub><mo stretchy="false">(</mo><mi>λ</mi><mn>,</mn><mi>u</mi><mn>,</mn><mi>ν</mi><mo stretchy="false">)</mo><mn>.</mn></mrow></math> Step 3. Hilbert filtering of g 2 ( λ , u , ν ) row by row <math display="block" id="M21"><mrow><msub><mi>g</mi><mn>3</mn></msub><mo stretchy="false">(</mo><mi>λ</mi><mn>,</mn><mi>u</mi><mn>,</mn><mi>ν</mi><mo stretchy="false">)</mo><mo>=</mo><mfrac><mn>1</mn><mrow><mn>2</mn><mi>π</mi></mrow></mfrac><msub><mi>g</mi><mn>2</mn></msub><mo stretchy="false">(</mo><mi>λ</mi><mn>,</mn><mi>u</mi><mn>,</mn><mi>ν</mi><mo stretchy="false">)</mo><mo>*</mo><msub><mi>h</mi><mi>H</mi></msub><mo stretchy="false">(</mo><mi>u</mi><mo stretchy="false">)</mo><mn>,</mn></mrow></math> where h H = 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 <math display="block" id="M22"><mrow><mi>f</mi><mo stretchy="false">(</mo><mover accent="true"><mi>r</mi><mo>→</mo></mover><mo stretchy="false">)</mo><mo>=</mo><mfrac><mn>1</mn><mn>2</mn></mfrac><mstyle displaystyle="true"><mrow><msubsup><mo>∫</mo><mrow><mi>θ</mi><mo>−</mo><mi>π</mi></mrow><mrow><mi>θ</mi><mo>+</mo><mi>π</mi></mrow></msubsup></mrow></mstyle><mfrac><mrow><mi>R</mi><mi>D</mi></mrow><mrow><msup><mrow><mrow><mo>[</mo><mrow><mi>R</mi><mo>−</mo><mover accent="true"><mi>r</mi><mo>→</mo></mover><mo>⋅</mo><msub><mover accent="true"><mi>e</mi><mo>→</mo></mover><mi>w</mi></msub><mo stretchy="false">(</mo><mi>λ</mi><mo stretchy="false">)</mo></mrow><mo>]</mo></mrow></mrow><mn>2</mn></msup></mrow></mfrac><msub><mi>g</mi><mn>3</mn></msub><mo stretchy="false">(</mo><mi>λ</mi><mn>,</mn><mi>u</mi><mn>,</mn><mi>ν</mi><mo stretchy="false">)</mo><mtext>d</mtext><mi>λ</mi><mn>.</mn></mrow></math> Mathematically, formula (4) is equivalent to formula (10) in reference [ 2 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 2 ( λ , 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 18 ]. When the weighted projection <math id="M23"><mrow><msubsup><mi>g</mi><mn>1</mn><mrow><mi>Z</mi><mi>o</mi><mi>u</mi></mrow></msubsup><mo stretchy="false">(</mo><mi>λ</mi><mn>,</mn><mi>u</mi><mn>,</mn><mi>ν</mi><mo stretchy="false">)</mo></mrow></math> is available, Zou et al . performed the Laplace operator based filtering by <math display="block" id="M24"><mrow><msubsup><mi>g</mi><mn>2</mn><mrow><mi>Z</mi><mi>o</mi><mi>u</mi></mrow></msubsup><mo stretchy="false">(</mo><mi>λ</mi><mn>,</mn><mi>u</mi><mn>,</mn><mi>ν</mi><mo stretchy="false">)</mo><mo>=</mo><mrow><mo>[</mo><mrow><mfrac><mrow><msup><mo>∂</mo><mn>2</mn></msup></mrow><mrow><mo>∂</mo><msup><mi>u</mi><mn>2</mn></msup></mrow></mfrac><mo>+</mo><mfrac><mrow><msup><mo>∂</mo><mn>2</mn></msup></mrow><mrow><mo>∂</mo><msup><mi>ν</mi><mn>2</mn></msup></mrow></mfrac></mrow><mo>]</mo></mrow><msubsup><mi>g</mi><mn>1</mn><mrow><mi>Z</mi><mi>o</mi><mi>u</mi></mrow></msubsup><mo stretchy="false">(</mo><mi>λ</mi><mn>,</mn><mi>u</mi><mn>,</mn><mi>ν</mi><mo stretchy="false">)</mo><mn>,</mn></mrow></math> and 2D Radon based filtering by <math display="block" id="M25"><mrow><msubsup><mi>g</mi><mi>F</mi><mrow><mi>Z</mi><mi>o</mi><mi>u</mi></mrow></msubsup><mo stretchy="false">(</mo><mi>λ</mi><mn>,</mn><mi>u</mi><mn>,</mn><mi>ν</mi><mo stretchy="false">)</mo><mo>=</mo><mo>−</mo><mfrac><mn>1</mn><mrow><mn>4</mn><msup><mi>π</mi><mn>2</mn></msup></mrow></mfrac><mstyle displaystyle="true"><mrow><msubsup><mo>∫</mo><mn>0</mn><mi>π</mi></msubsup></mrow></mstyle><mo>|</mo><mi>cos</mi><mi>μ</mi><mo>|</mo><msubsup><mi>g</mi><mn>3</mn><mrow><mi>Z</mi><mi>o</mi><mi>u</mi></mrow></msubsup><mo stretchy="false">(</mo><mi>λ</mi><mn>,</mn><mi>u</mi><mn>,</mn><mi>ν</mi><mo stretchy="false">)</mo><mtext>d</mtext><mi>μ</mi><mn>,</mn></mrow></math> where <math display="block" id="M26"><mrow><msubsup><mi>g</mi><mn>3</mn><mrow><mi>Z</mi><mi>o</mi><mi>u</mi></mrow></msubsup><mo stretchy="false">(</mo><mi>λ</mi><mn>,</mn><msup><mi>u</mi><mo>*</mo></msup><mn>,</mn><msup><mi>ν</mi><mo>*</mo></msup><mo stretchy="false">)</mo><mo>=</mo><mstyle displaystyle="true"><mrow><msub><mo>∫</mo><mrow><msub><mi>Ω</mi><mi>λ</mi></msub></mrow></msub></mrow></mstyle><msubsup><mi>g</mi><mn>2</mn><mrow><mi>Z</mi><mi>o</mi><mi>u</mi></mrow></msubsup><mo stretchy="false">(</mo><mi>λ</mi><mn>,</mn><mi>u</mi><mn>,</mn><mi>ν</mi><mo stretchy="false">)</mo><mi>δ</mi><mo stretchy="false">(</mo><mi>u</mi><mi>cos</mi><mi>μ</mi><mo>+</mo><mi>ν</mi><mi>sin</mi><mi>μ</mi><mo>−</mo><mi>s</mi><mo stretchy="false">)</mo><mtext>d</mtext><mi>u</mi><mtext>d</mtext><mi>μ</mi><mn>,</mn></mrow></math> 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 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. Also, the truncated 2D Shepp-Logan phantom is filtered by the three kinds of filters and the results are shown in Fig. 3 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. 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 1 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 2 In the following experiment, we demonstrated the efficiency of our method on simulated projection data from the 3D Shepp-Logan phantom [ 20 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 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 TABLE 1 . Because the width of the detector is 256 pixels in the Simulation 2, the simulated projections are laterally truncated. 1 1 parameters 1 1 Simulation 1 1 1 Simulation 2 1 1 source to 1 1 150 mm 1 1 150 mm 1 1 object distance 1 1 1 1 1 1 source to 1 1 600 mm 1 1 600 mm 1 1 detector distance 1 1 1 1 1 1 size of detector 1 1 256 pixel×256 pixel 1 1 256 pixel×512 pixel 1 1 size of detector cell 1 1 0.2 mm×0.2 mm 1 1 0.2 mm×0.2 mm 1 1 number of 1 1 360 1 1 360 2 1 projections per turn 1 1 1 1 1 1 pitch of helix 1 1 25 mm 1 1 25 mm 1 1 Number of turn 1 1 2 1 1 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. 4 . Fig. 5 Fig. 5 illustrates the profiles of middle row data of the first row in Fig. 4 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 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. 1 1 Simulation 1 1 The helical FDK algorithm 1 1 LORA method 1 1 The proposed method 1 1 Simulation 1 1 1 64 1 1 9232 1 1 92 1 1 Simulation 2 1 1 38 1 1 5666 1 1 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 15 ]. The reconstruction images from the simulated truncated projections are shown in Fig. 6 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 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 TABLE 2 . The conclusions are similar to Simulation 1. B. Real data 2 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 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 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 Figs. 8 and 9 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 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. V. CONCLUSION 1 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 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 http://doi.org/10.3233/XST-130406 Practical geometric calibration for helical cone-beam industrial computed tomography G Wang , T H Lin , P C Cheng , et al . A general cone-beam reconstruction algorithm . IEEE T Med Imaging , 1993 , 12 : 486 - 496 . DOI: 10.1109/42.241876 http://doi.org/10.1109/42.241876 A general cone-beam reconstruction algorithm A Katsevich . A general scheme for constructing inversion algorithms for cone beam CT . Int J Math Math Sci , 2003 , 21 : 1305 - 1321 . DOI: 10.1155/S0161171203209315 http://doi.org/10.1155/S0161171203209315 A general scheme for constructing inversion algorithms for cone beam CT Y Zou , X Pan and E Y Sidky . Theory and algorithms for image reconstruction on chords and within regions of interest . J Opt Soc Am A , 2005 , 22 : 2372 - 2384 . DOI: 10.1364/JOSAA.22.002372 http://doi.org/10.1364/JOSAA.22.002372 Theory and algorithms for image reconstruction on chords and within regions of interest F Natterer . The Mathematics of Computerized Tomography . American Society for Industrial & Applied Mathematics Press , 2001 , 121 - 125 . DOI: 10.1137/1.9780898719284 http://doi.org/10.1137/1.9780898719284 The Mathematics of Computerized Tomography K S Sharma , C Holzner , D M Vasilescu , et al . Scout-view assisted interior micro-CT . Phys Med Biol , 2013 , 58 : 4297 - 4314 . DOI: 10.1088/0031-9155/58/12/4297 http://doi.org/10.1088/0031-9155/58/12/4297 Scout-view assisted interior micro-CT H Y Yu and G Wang . Compressed sensing based interior tomography . Phys Med Biol , 2009 , 54 : 2791 - 2805 . DOI: 10.1155/2009/125871 http://doi.org/10.1155/2009/125871 Compressed sensing based interior tomography J Yang , H Yu , M Jiang , et al . High order total variation minimization for interior tomography . Inverse Probl , 2010 , 26 : 035013 . DOI: 10.1088/0266-5611/26/3/035013 http://doi.org/10.1088/0266-5611/26/3/035013 High order total variation minimization for interior tomography X Jin , A Katsevich , H Yu , et al . Interior tomography with continuous singular value decomposition . IEEE T Med Imaging , 2012 , 31 : 2108 - 2119 . DOI: 10.1109/TMI.2012.2213304 http://doi.org/10.1109/TMI.2012.2213304 Interior tomography with continuous singular value decomposition F Dennerlein . Cone-beam ROI reconstruction using the Laplace operator . 11th international meeting on Fully Three-dimensional image reconstruction in Radiology and Nuclear Medicine , 2011 , 80 - 83 . Cone-beam ROI reconstruction using the Laplace operator Y Xia , A Maier , F Dennerlein , et al . Efficient 2D filtering for cone-beam VOI reconstruction , Nuclear Science Symposium and Medical Imaging Conference (NSS/MIC) , Anaheim, CA , October 27, 2012 , 2415 - 2420 . DOI: 10.1109/NSSMIC.2012.6551549 http://doi.org/10.1109/NSSMIC.2012.6551549 Efficient 2D filtering for cone-beam VOI reconstruction F Dennerlein and A Maier . Region-of-interest reconstruction on medical C-arms with the ATRACT algorithm , SPIE Medical Imaging 2012: Physics of Medical Imaging , San Diego, CA , February 2012 , 83131B . DOI: 10.1117/12.913274 http://doi.org/10.1117/12.913274 Region-of-interest reconstruction on medical C-arms with the ATRACT algorithm Y Xia , A Maier , H G Hofmann , et al . Reconstruction from truncated projections in cone-beam CT using an efficient 1D filtering , SPIE Medical Imaging 2013: Physics of Medical Imaging , Lake Buena Vista, FL , February 2013 , 86681C . DOI: 10.1117/12.2007484 http://doi.org/10.1117/12.2007484 Reconstruction from truncated projections in cone-beam CT using an efficient 1D filtering F Dennerlein and A Maier . Approximate truncation robust computed tomography-ATRACT . Phys Med Biol , 2013 , 58 : 6133 - 6148 . DOI: 10.1088/0031-9155/58/17/6133 http://doi.org/10.1088/0031-9155/58/17/6133 Approximate truncation robust computed tomography-ATRACT Y Xia , H Hofmann , F Dennerlein , et al . Towards clinical application of a Laplace Operator-based region of interest reconstruction algorithm in C-Arm CT . IEEE T Med Imaging , 2014 , 33 : 593 - 606 . DOI: 10.1109/TMI.2013.2291622 http://doi.org/10.1109/TMI.2013.2291622 Towards clinical application of a Laplace Operator-based region of interest reconstruction algorithm in C-Arm CT L A Feldkamp , L C Davis and J W Kress . Practical cone-beam algorithm . J Opt Soc Am A , 1984 , 1 : 612 - 619 . DOI: 10.1364/JOSAA.1.000612 http://doi.org/10.1364/JOSAA.1.000612 Practical cone-beam algorithm M Grass , T Köhler and R Proksa . 3D cone-beam CT reconstruction for circular trajectories . Phys Med Biol , 2000 , 45 : 329 - 347 . DOI: 10.1088/0031-9155/45/2/306 http://doi.org/10.1088/0031-9155/45/2/306 3D cone-beam CT reconstruction for circular trajectories X B Zou , H Yu and L Zeng . Laplace operator based reconstruction algorithm for truncated spiral cone beam computed tomography . J X-ray Sci Technol , 2013 , 21 : 515 - 526 . DOI: 10.3233/XST-130398 http://doi.org/10.3233/XST-130398 Laplace operator based reconstruction algorithm for truncated spiral cone beam computed tomography X C Wang , B Yan , L Li , et al . Cone-beam local reconstruction based on a Radon inversion transformation . Chinese Phys B , 2012 , 21 : 8702 . DOI: 10.1088/1674-1056/21/11/118702 http://doi.org/10.1088/1674-1056/21/11/118702 Cone-beam local reconstruction based on a Radon inversion transformation P A Toft . The Radon transform, theory and implementation. Ph.D. Thesis , Technical University of Denmark , 1996 . The Radon transform, theory and implementation. Ph.D. Thesis 10.13538/j.1001-8042/nst.26.020401 Supported by the National High Technology Research and Development Program of China (No. 2012AA011603) and National Nature Science Foundation of China (No. 61372172) 2014-04-11 2014-05-15 2015-04-20 2015-04-20 2021-04-11 2015 Nuclear Science and Techniques 2 26