Three-dimensional weighting reconstruction algorithm for circular cone-beam CT under large scan angles

LOW ENERGY ACCELERATOR, RAY AND APPLICATIONS

Three-dimensional weighting reconstruction algorithm for circular cone-beam CT under large scan angles

Ya-Fei Yang，

Ding-Hua Zhang，

Kui-Dong Huang，

Fu-Qiang Yang，

Zong-Zhao Gao

Nuclear Science and Techniques

Vol.28, No.8

Article number 116

Published in print 01 Aug 2017

Available online 07 Jul 2017

DOI：10.1007/s41365-017-0262-3

55203

Improving imaging quality of cone-beam CT under large cone angle scan has been an important area of CT imaging research. Considering the idea of conjugate rays and making up missing data, we propose a three-dimensional (3D) weighting reconstruction algorithm for cone-beam CT. The 3D weighting function is added in the back-projection process to reduce the axial density drop and improve the accuracy of FDK algorithm. Having a simple structure, the algorithm can be implemented easily without rebinning the native cone-beam data into cone- parallel beam data. Performance of the algorithm is evaluated using two computer simulations and a real industrial component, and the results show that the algorithm achieves better performance in reduction of axial intensity drop artifacts and has a wide range of application.

FDK algorithmMissing dataConjugate raysThree-dimensional weighting

1. Introduction

With notable features of high efficiency, high precision, smaller size, smaller radiation dose and lower cost, cone-beam computed tomography (CBCT) has been widely used in medical imaging and industrial non-destructive testing [1, 2]. Currently in CBCT systems of circular scanning geometry, the FDK (Feldkamp–Davis–Kress) algorithm, a straightforward generalization of the 2D fan-beam geometry to 3D cone-beam geometry, is the most popular method [3]. As an analytical algorithm, FDK algorithm is advantageous in its simple structure, low calculation expense, easy implementation and accuracy [4, 5]. However, as an inexact algorithm, FDK algorithm may introduce significant artifacts in the images [6]. According to the data sufficiency condition (DSC), the circular scanning geometry, due to its own structure, always has the problem of a donut-shaped area of missing data in Radon space [7], and because of the data-insufficiency, FDK algorithm fails to reconstruct the image accurately [8]. An increase of cone angle in FDK algorithm can cause two kinds of problems: high-frequency cone-beam artifacts and axial intensity drop. The former makes the images deviate from the real object greatly. Iterative algorithm may solve this image quality issues at wider cone angles, but it cannot meet the actual needs since large-scale matrix operations require a large memory and a particularly time-consuming process when reconstructing high-resolution images (such as 2048³ or 4096³, which is popular in industrial CT).

Many methods have been proposed to reduce cone-beam artifacts. They can be classified into four categories:

1) Satisfying the DSC condition by turning the circular trajectory into other scanning trajectories, such as ‘circle + circle’ [9],　‘circle + line’ [10],　‘circle + spiral’ [11];

2) Modifying FDK by rebinning the native cone-beam data into cone-parallel beam data, such as T-FDK [12], HT-FDK [13], CW-FDK [14], CB-FBP [15], ACE [16, 17], C-FDK [18],　BPW-FDK [19];

3) Adding a correction term in the FDK algorithm to get the information contained in a circular cone-beam scan but not utilized in the FDK algorithm based on Grangeat formula [20]; and

4) Modifying FDK by adding a weighting function in the reconstruction process, such as xFDK [21], FAFDK [22], Weighted FDK [23, 24].

However, the Category 1 methods require a high precision trajectory and a large amount of computation, with increased difficulty in mechanical design, hence the low feasibility in practical applications. The Category 2 methods, in which the cone-beam has to be rebinned into cone-parallel geometry, may cause spatial resolution loss especially when cone-beam data samples are sparse, and they are useful just when the cone angle is small (≤±5°) [12-19, 25]. The Category 3 methods are more efficient with less intensity drop than FDK, but need huge amount of calculation [26-28]. The Category 4 methods have simple and intuitive structures, with reduced cone-beam artifacts (especially for intensity drop in the z axis) [21-24].

This paper presents a new algorithm based on FDK. Based on the Tang’s idea to make up missing data in Radon space, a 3D back-projection weighting function is introduced into FDK to reduce the reconstruction artifacts [15].　Two phantoms and an industrial component are used to test feasibility and effectiveness of the algorithm.

2. Materials and methods

2.1 Cone-beam scanning and missing data

The circular source trajectory of FDK is shown in Fig. 1 (we mainly discuss the flat-panel detector in this paper), where the coordinate system is O-xyz; S is the X-ray source; D is the virtual flat-panel detector; R is radius of the circular source trajectory; and P(x,y,z) is a voxel on reconstruction objective f(x,y,z) which can be determined uniquely by the rotation angle θ, fan angle γ and cone angle α. The FDK algorithm can be expressed as:

Fig. 1.

(Color online) Schematic diagram of the native CB geometry.

f (x, y, z) = \int_{0}^{2 π} \frac{d s o^{2}}{{(R + x \cos θ + y \sin θ)}^{2}} \times \frac{1}{2} \times (\frac{d s o}{\sqrt{R^{2} + u^{2} + v^{2}}} g (θ; u, v)) \otimes h (u) d θ,

(1)

where, g(θ,u,v) is the projection data obtained from D, (u,v) is coordinate value of the detector system, $\otimes$ is the row-wise 1D convolution and h(u) is the ramp filtering kernel [3].

FDK algorithm achieves good results only under small cone angle scan, but sometimes one has to reconstruct CT images under large cone angle scan. For example, the source-detector distance of micro-CT is much smaller than that of normal CT (so as to obtain enough signals from the detector of low power output), and this means a wider cone angle scan. In this case, the cone-beam artifacts especially for the intensity drop become even more significant with the increase of cone angle.

An approach to reduce high-frequency cone-beam artifacts is the idea of conjugate rays (For any ray passing through the voxel of P, there only exists a corresponding conjugate ray, and the rotation angle between the direct ray and the conjugate ray is always 180°). The inconsistency of conjugate rays causes the cone-beam artifacts for FDK algorithm, which treats each pair of conjugate rays equally with a fixed weight of 0.5 [15]. To reduce cone-beam artifacts, the ray being farther away from the source S and passing through the voxel P(x,y,z) shall have greater weight in the projection in a pair of conjugate rays (the smaller the cone angle, the more reliable the projection), and the other one has a smaller weight.

Another approach to reduce cone-beam artifacts is to make up missing data (mainly for the axial intensity drop). According to the DSC condition, as shown in Fig. 2, the Radon space data (the first derivative of the Radon transform) are not sufficient in the z direction for the circular trajectory [28]. With the increase of z, the shaded area increases gradually, forming a structure similar to a donut [7].　However, FDK algorithm fills the shaded area with exactly zeroes, which results in missing information.

Fig. 2.

Radon missing data of circular scanning.

For the problem of missing data in FDK, studies show that deviation in the z direction can be described as a hilltop-like function (in which the first derivative is an increasing function)　and it is feasible to compensate missing data in the z direction using typical experiential functions, such as cosine function and Gaussian function [23]. At this time the weight in the back-projection process is equivalent to an experiential hilltop-like function, rather than a fixed weight of 0.5 in FDK.

2.2 3D weighted reconstruction under large cone angle

The algorithm in Section 2.1 can be described from different aspects especially for the second and fourth methods. Here we take CB-FBP and Weighted FDK as examples. From the idea of conjugate ray, the experiential 3D weight of CB-FBP is：

w_{3 d}^{C B - F B P} (α, z) = \frac{1}{1 + {(\tan (| α |) / \tan (| α_{c} |))}^{k | z |}},

(2)

and it meets the normalization condition:

w_{3 d}^{C B - F B P} (α, z) + w_{3 d}^{C B - F B P} (α_{c}, z) = 1 .0,

(3)

where, α and α_c are the cone angle corresponding to a direct ray and the conjugate ray, respectively; and |k| is an experiential function increasing with |z|. Starting from making up missing data, the experiential 3D weight of Weighted FDK is:

w_{3 d}^{W e i g h t e d F D K} (z, c_{1}, c_{2}) = \frac{1}{2 \times \cos [c_{1} z / (R - c_{2} r)]},

(4)

and it satisfies the inequality condition:

w_{3 d}^{W e i g h t e d F D K} (z, c_{1}, c_{2}) + w_{3 d}^{W e i g h t e d F D K} (z, c_{1}, c_{2}) \geq 1 .0,

(5)

where, r=(x2+y+2z2)1/2, c₁ and c₂ are experiential parameters provided for accelerating or decelerating the z distance variable and the radial distance variable r, respectively, thereby adjusting the shape of the weight function [23]. For the hilltop-like description, the selection of parameters c₁ and c₂ should satisfy mathematical constraint conditions:

0 \leq c_{1} z / (R - c_{2} r) \leq π / 2,

(6)

However, mathematically, there is a big weakness in Weighted FDK: There is no doubt than we can decrease the axial intensity drop at z1 for normal detector once we have chosen c₁ and c₂; however, there must be a z₂ for larger detector which ruins the mathematical constraint conditions, leading to

$\lim_{z \to z_{2}^{-}} w_{3 d}^{W e i g h t e d F D K} (z) \to + \infty or \lim_{z \to z_{2}^{+}} w_{3 d}^{W e i g h t e d F D K} (z) \to - \infty .$

This means the gray values may be +∞ or −∞ at z₂ for larger detector even if the gray value is close to zero (the gray value will never be zero in simulations or practice), will cause great distortion in the image. This forces us to re-select the c₁ and c₂ when we scan the same object using the same scanning parameters without any other change (We just turn a normal detector into a larger detector). The application of Weighted FDK is greatly restricted.

Both of CB-FBP and Weighted FDK can be expressed in a general formula as follows:

f (x, y, z) = \int_{0}^{2 π} \frac{d s o^{2}}{{(R + x \cos θ + y \sin θ)}^{2}} \times w_{3 d} \times (\frac{d s o}{\sqrt{R^{2} + u^{2} + v^{2}}} g (θ; u, v)) \otimes h (u) d θ,

(7)

where w_3d is a 3D weighting function, and the method to determine w_3d differs from algorithm to algorithm, and so do the starting idea.

Mathematically, we can analyze the differences between the two algorithms and FDK. For FDK, the two weights’ sum of conjugate rays is always 1.0, and the two weights are always 0.5; while for CB-FBP, the two weights are always not the same in the non-central plane, though the two weights’ sum of conjugate rays is always 1.0. On the contrary, for Weighted FDK, the two weights of conjugate rays are always the same; but their sum is always bigger than 1.0 in the non-central plane. The differences in mathematics may be their own advantages over FDK. Since CB-FBP is designed against the high-frequency cone-beam artifacts without any consideration of the axial intensity drop (which is more significant), in contrast, Weighted FDK (which introduces an ad hoc correction) is more effective for axial intensity drop. The inequality condition in Weighted FDK is believed to be more useful for axial intensity drop than the normalization condition in CB-FBP. However, CB-FBP is useful in small cone angle only and may cause spatial resolution loss. Besides, Weighted FDK does not have a widespread compatibility and may not be applicable in practice especially for complex objects since it has two parameters to adjust through experience and a big weakness in mathematics [18, 29].

In this paper, we have a comprehensive consideration of CB-FBP and Weighted FDK. Starting from CB-FBP (conjugate rays) and Weighted FDK (making up missing data), the cone-beam differs from fan-beam in the cone angle itself (the cone angle may cause cone-beam artifacts) and a 3D weighting function contacted with cone angle is added in the back- projection process to decrease cone-beam artifacts. Obviously, the modifying function should meet six basic conditions:

1) The modifying weight should be related to cone angle, and symmetrical about the central plane;

2) The bigger the cone angle, the bigger the modifying function value;

3) The new proposed algorithm should be equivalent to FDK at the central plane (the cone angle is zero);

4)The weighting function should satisfy the inequality condition of Weighted FDK in order to make up missing data or reduce axial intensity drop;

5) The weighting function should be a hilltop-like function just like Weighted FDK; and

6)The new proposed algorithm does not have big weakness in mathematics.

According to the geometry in Fig.3, the loss of information has much to do with the increase of the cone angle in the back-projection process, and the loss rate is determined by the cone angle from the source S to the voxel P₁(x,y,z). We can add a weighting function in the back-projection progress to make up the loss of information.

Fig. 3.

The cone angle under the native CB geometry.

Therefore, we define the 3D modifying weight as:

w_{3 d} (z, t, s) = \frac{1}{2} \sqrt{1 + p * \tan^{2} (α)} = \frac{1}{2} \sqrt{1 + \frac{p * z^{2}}{{(R + s)}^{2} + t^{2}}}, p \geq 0,

(8)

where, t = ycosθ – xsinθ, s = xcosθ + ysinθ and p is an experiential parameter to be adjusted like other algorithms (the proposed algorithm is the same as FDK when p =0, and the number of p means the degree of compensation in the z direction). Eq.(8) can be rewritten as:

w_{3 d} (x, y, z) = \frac{1}{2} \sqrt{1 + \frac{p * z^{2}}{{(R + x \cos θ + y \sin θ)}^{2} + {(y \cos θ - x \sin θ)}^{2}}},

(9)

Different from FDK, in the proposed algorithm, the weights of conjugate rays are always unequal and bigger than 0.5 in the non-central plane. It does not have the big weakness of Weighted FDK (the weight will never be +∞ or −∞), and it is a hilltop-like function and meets the fourth condition:

w_{3 d} (z, t, s) + w_{3 d} (z, t, s_{c}) \geq 1 .0,

(10)

The degree of compensation in the z direction can be expressed as:

w_{3 D}^{m a k i n g u p} = w_{3 d} (z, t, s) + w_{3 d} (z, t, s_{c}) -1 = \frac{1}{2} \sqrt{1 + \frac{p * z^{2}}{{(R + s)}^{2} + t^{2}}} + \frac{1}{2} \sqrt{1 + \frac{p * z^{2}}{{(R - s)}^{2} + t^{2}}} -1, p \geq 0,

(11)

Figure 4 shows the 3D weight as functions of t and s, at p =1.0, R =375 mm, and z =±100 mm. The surface shape is similar to a half saddle, proving that the function is also a hilltop-like function. The surface is relatively high in the middle and relatively low on both sides in the t direction, and reduces as s increases (For the voxel in a plane, the farther away from the source, the smaller the cone angle, and the smaller the weight).

Fig. 4.

(Color online) The 3D weight as functions of t and s.

Finally, the proposed algorithm starting from making up missing data and conjugate rays can be expressed as:

f (x, y, z) = \int_{0}^{2 π} \frac{d s o^{2}}{{(R + x \cos θ + y \sin θ)}^{2}} \times w_{3 d} (x, y, z) \times (\frac{d s o}{\sqrt{R^{2} + u^{2} + v^{2}}} g (θ; u, v)) \otimes h (u) d θ,

(12)

3. Results and Discussion

3.1 The experiment object and parameters

Simulations were performed with the standard 3D Shepp-Logan phantom [30] and Defrise disk phantom [31] to evaluate FDK, Weighted FDK and the proposed algorithm. An industrial component (the data were obtained from the authors’ laboratory) was tested to evaluate the three algorithms. The experiment parameters are given in Table 1.

System parameters of CBCT in simulations and experiment (scan angle,360°; detector pixels, 512×512; image matrix, 512×512×512).

Parameters	3D Shepp-Logan	Defrise disk	Industrial object
Source-to-rotation-center distance (mm)	480	480	978.9
Source-to-detector-center distance (mm)	960	960	1235.7
Detector size (mm)	512×512	512×512	102.4×102.4
Size of each voxel (mm)	0.5	0.5	0.16
Cone-angle of the X-ray beam (°)	±15	±15°	±2.4°

In order to get correct parameters in Weighted FDK, we chose c₂=0 first and adjusted c₁ to correct for the intensity drop along x= y= z. This determined c₁ and tuned c₂ to correct for residual drop as a function of radial distance from the origin. Besides, we added the image of fan-beam CT scans of the same object as a reference in two simulations.

3.2 Reconstruction results

Figure 5 shows reconstruction results of the standard 3D Shepp-Logan phantom at y=−25 mm, c₁=1.32, c₂=0.05 and p=1.87. While the FDK results were the worst, both Weighted FDK and the proposed algorithm achieved good results. However, adjusting the two parameters in Weighted FDK was more complicated than in the proposed algorithm.

Fig. 5.

(Color online) Reconstructions of the standard 3D Shepp-Logan phantom, x-z views. Display window [0.98, 1.05]: (a) fan-beam CT scan, (b) FDK, (c) Weighted FDK, (d) Proposed algorithm, (e) The centrally vertical profiles of the slices in Figs. 5(b), 5(c) and 5(d).

Figure. 6 shows reconstruction results of the Defrise disk phantom at y=0 mm, c₁=4.8, c₂=0.2 and p=120. The FDK did badly. The Weighted FDK improved the image quality in some degree, but the intensity drop still existed in the z direction (there would be a distortion in the image if we enlarged c₁ or c₂). With a better mathematical model, the proposed algorithm eliminated the intensity drop to a large degree, and made great corrections to the subjective sight and the profiles.

Fig. 6.

(Color online) Reconstructions of the Defrise disk phantom, x-z views. Display window [0.5, 1.1]: (a) fan-beam CT scan, (b) FDK, (c) Weighted FDK, (d) Proposed algorithm, (e) The centrally vertical profiles of slices in Figs. 6(b), 6(c) and (d).

Figure 7 shows reconstruction results of the industrial component at y=0 mm, c₁=32, c₂=0.3 and p=3000 (There was a small image distortion by enlarging c₁ or c₂). Comparing with FDK, the correction of Weighted FDK algorithm was inconspicuous, while the proposed algorithm performed the best, with good subjective sight and profiles.

Fig. 7.

(Color online) Reconstructions of the industrial component, x-z views. Display window [0.5, 1.5]: (a) FDK, (b) Weighted FDK, (c) Proposed algorithm, (d) The profiles by the three algorithms on the location marked as line C in Fig. 7(a).

Table 2 is the numerical comparison of slice quality in Fig. 7, where rectangle A is the comparison position of the signal to noise ratio (SNR), and rectangle B is the comparison position of the contrast to noise ratio (CNR) and average gradient (AG) [32]. It can be seen clearly from Table 2 that SNR and AG of the proposed algorithm are much better than that of the Weighted FDK, but the improvement in CNR by this algorithm is the same as Weighted FDK because of increases in stochastic, quantum noise of the reconstruction as a result of different weights assigned to the rays (A noise-optimal reconstruction algorithm like standard FDK would assign equal weight to the two rays).

Numerical comparison of image quality (the bracketed percentages are improvements to the FDK values)

Algorithm	SNR	CNR	AG
FDK	5.596	3.281	0.0404
Weighted FDK	6.688(19.5%)	3.296(0.46%)	0.0463(14.6%)
This algorithm	7.340(31.2%)	3.295(0.43%)	0.0534(32.2%)

From the results of the standard 3D Shepp-Logan phantom, Defrise disk phantom and the industrial component, it is easy to find FDK, Weighted FDK and the proposed algorithm are completely equivalent at the center plane. Weighted FDK is effective when the cone beam artifacts are inconspicuous; but it cannot work well when the artifacts became conspicuous, as the mathematical model may be not suitable. With a better mathematical model Combining the CB-FBP algorithm and Weighted FDK algorithm, the proposed algorithm achieves better results in the simulations and the real objects and is more valuable in practice, because the artifacts of industrial objects are more conspicuous than simulations.

4. Conclusion

From the view of conjugate rays and making up missing data, we present a 3D weighting cone-beam CT reconstruction algorithm under large cone angle scan. Keeping in circular source trajectory, without any action of rebinning, the algorithm just adds a weighting function in the back-projection progress and can be realized easily, with the data flow being identical to FDK. The algorithm is equivalent to FDK in the central plane. Improving the accuracy of original FDK significantly and reducing the cone angle artifacts, the proposed algorithm achieves better image quality than FDK and Weighted FDK in phantom simulation and real object test. Besides, it has good adaptability under large cone angle scan (±15°). For real objects, the algorithm improves image quality greatly even in small cone angle (because artifacts of real objects are more conspicuous than simulation phantoms), indicating its value in practical applications. The parameter-adjustment time can be saved as one needs to adjust just one parameter, p. For super complex objects, p can be an experiential function if necessary. The algorithm can be widely applied in medical diagnosis and non-destructive testing.

References

J. Hao, L. Zhang, L. Li et al.,

A practical image reconstruction and processing method for symmetrically off-center detector CBCT system

. Nucl. Sci. Tech. 24(4), 40204-040204 (2013). doi: 10.13538/j.1001-8042/nst.2013.04.004