# Reconstruction Method for Nonconcyclic Dual-Source Circular Cone-Beam CT with a Large Field of View

Authors: Ming Chen; Gang Li; Weiwei Qi; Jing Zou; Yong-guo Zheng
Journal: Mathematical Problems in Engineering
 (2013)

---

## Abstract

In industrial computed tomography (CT), it is often required to inspect large objects whose size is beyond a reconstructed field of view (FOV). Some multiscan modes have been developed to acquire the complete CT projection data for a larger object using small panel detectors. In this paper, we give a non-concyclic dual-source circular cone-beam scanning geometry based on the idea of multiscan modes and propose a backprojection-filtration-based (BPF) reconstruction algorithm without data rebinning. Since the FOV calculated according to this nonconcyclic dual-source circular CT scanning geometry is larger than cardiac dual-source CT scanning geometry, our method can reconstruct larger horizontal slices (i.e., the slices perpendicular to rotation axis) than cardiac dual-source CT. The quality of CT images is expected to be superior to those obtained using larger panel detectors. The simulation results have indicated that CT images obtained by the proposed method are satisfying.


---

## Body

## 1. Introduction

In industrial CT, it is often required to inspect large objects whose size is beyond a field of view (FOV) calculated according to CT scanning geometry, where the maximum horizontal (i.e., perpendicular to rotation axis) FOV is usually fixed in CT system. When a slice of the inspected object cannot be completely covered within the FOV of the cone beam, complete projection data passing through the object cannot be obtained. However, conventional CT reconstruction algorithm, for example, filtered backprojection (FBP), requires the complete projection data. Therefore, those conventional scan modes and reconstruction methods can not be used. To solve this problem, several multiscan modes have been developed to acquire the complete projection data. Of them, the rotation-translation-translation (RTT) multiscan mode [1, 2] and rotation-translation (RT) multiscan mode [3, 4] are typical. In fact, RT and RTT multiscan modes not only acquire completed projection data, but can also reduce the differences of the flux intensities detected by different detector cells.For fan-beam RT multiscan modes, the reconstruction algorithm based on backprojection- filtration (BPF) was proposed by Chen et al. [5, 6], which can exactly reconstruct the images and does not introduce data rebinning. BPF algorithm, as the basis of our work, was first proposed by Zou and Pan [7] to exactly reconstruct the images for helical cone beam CT. After that, an explicit BPF formula for 2D image reconstruction was proposed by Noo et al. [8] using the finite inversion formula of Hilbert transform.In this work, we give a nonconcyclic dual-source circular cone beam scanning geometry based on the idea of RT multiscan scan mode and propose a BPF-based reconstruction formula. Since The FOV calculated according to this nonconcyclic dual-source circular CT scanning geometry is larger than cardiac dual-source CT scanning geometry, our method can reconstruct larger horizontal slices. The quality of CT images reconstructed is expected to be superior to those obtained using larger panel detectors. The simulation results have indicated that CT images obtained by the proposed algorithm are satisfying. Since we make use of the approximate idea of Feldkamp (FDK) algorithm [9] in reconstructing formula deduction, reconstruction images are satisfying for a small cone angle, which are usually better for static objects. For dynamic process, there are some results in [10, 11].


## 2. Nonconcyclic Dual-Source Circular Scan Mode with a Larger FOV

In this section, we give a nonconcyclic dual-source circular cone beam scan mode with a larger FOV based on RT dual-scan modes.As shown in Figure1, there are two pairs of X-ray sources and panel detectors. Let S1 and S2 denote the focuses of X-ray sources. Let D1 and D2 denote panel detectors. In the nonconcyclic dual-source circular scan mode, two pairs of X-ray sources and panel detectors are fixed, and the inspected object is placed on the turntable. Now, we give this nonconcyclic dual-source scanning geometry. Firstly, the center of the turntable O is located in the cone beam formed by S1 and D1. For some large inspected objects, one cone beam formed by S1 and D1 cannot completely cover its horizontal slice. So, secondly, we set S2 and D2 making S1, S2 and D1, D2 located at the same side of O, and two pairs of cone beams are parallel, as shown in Figure 1. Now, we give some relevant parameters to describe this scanning geometry. Let R2 denote the distance from S1 to D1 (equal to the distance from S2 to D2). Let R1 denote the distance from S1 to the virtual panel detectors which is through O and parallel to panel detectors D1. Two pairs of panel detectors have the same size, where L and W denote the length and width of the panel detectors, respectively. To obtain sufficient projection data, the distance from S2 to S1 must be less than LR1/(2R2) [6]. According to the location of the three points O, S1, and S2, we know that distance from S2 to O is larger than the distance from S1 to O. So two X-ray sources are nonconcyclic in this scanning geometry. The maximum horizontal FOV is 4L(R2-D)/(4L+4R22+L2) in nonconcyclic scanning geometry [6], where D is the distance from collimator of the detector to the crystal of the detector.Figure 1
The nonconcyclic dual-source circular cone beam scanning geometry with a larger FOV.


## 3. Reconstruction Formula

For the nonconcyclic dual-source circular scan mode above, we give the reconstruction formula in this section. We know that the key point to BPF algorithm is to obtain the Hilbert image from projections along some directions. The problem we face is how to weigh each set of the derivatives of two pairs of projection data to obtain two differentiated backprojection (DBP) images, and how to merge them into an entire DBP image that is related to the Hilbert image.For the deduction of the reconstruction formula, we need give some denotations for the scan mode in Section2. As shown in Figure 2, we establish 3D Cartesian coordinate system Oξηz, where z-axis is the rotating axis. Like the deduction of many CT reconstruction formulas, we also adopt two virtual panel detectors whose centers are O1 and O2, where the three points O1, O2, and O are collinear. Let O1′,O2′ denote the projection point of S1,S2 to the upper side line of two virtual panel detectors, as shown in Figure 2. We need to establish rotating coordinate system Oxyz in the deduction, and let β denote the rotating angle that is a clockwise angle from ξ-axis to x-axis. We set the 2D coordinate system Oiuivi (i=1,2) in two virtual panel detectors, where the direction vector of ui in Oxyz coordinate system is β=(cosβ,sinβ,0), and the direction vector of vi is the same direction with z-axis.Figure 2
Definition of coordinate systems and geometrical parameters.Letf(x0) denote the density function of the inspected object where x0=(x0,y0,z0) is a reconstructed point on the inspected object. Let bθ(x0) denote the DBP image of f(x0), where θ is an angle measured from the y-axis anticlockwise in the plane Oxy. Let pi(β,ui,vi) denote the projection data under X-ray source Si, and bθ,i(x0) denote the DBP image from pi(β,ui,vi). Let hi denote the distance between the points O and Oi, and hi must satisfy 0<h1<R1L/(2R2) and h2<3R1L/(2R2) [6] that indicates that the projection area of the inspected objects under two X-ray sources is some what overlap. Let (u0,iv0,i) denote the projection address of x0 in Oiuivi system under X-ray source Si. Let i=1,2 above.Now, using the idea of FDK algorithm, we deduct the formula of DBP imagebθ(x0) from the fan-beam DBP formula for RT multiscan modes [5, 6]. The formula deduction is divided into two kinds of situations: (i) x0=(x0,y0,z0) is on the middle plane where z0=0; (ii) x0=(x0,y0,z0) is on the off-middle plane where z0≠0.(i) Ifx0 is on the middle plane, we may directly obtain the cone beam DBP image bθ(x0) from the fan-beam DBP [6] as follows:
(1)bθ(x0)=bθ,1(x0)+bθ,2(x0),
where
(2)bθ,1(x0)=∫02πR1(R1-x0·β⊥)2ddu1×{[kε(R1(u1+h1)R12+u12)-kε(R1(u1+h1)R12+u12-T1)]×R12-h1u1R12+u12×p1(β,u1,0)×sgn(×u1R1-θsin(×u1R1-θβ+tan-1×sgnsinsin×u1R1-θ))}|u1=u0,1dβ,(3)bθ,2(x0)=∫02πR1(R1-x0·β⊥)2ddu2×{(kε(R1(u2+h2)R12+u22-T1))×R12-h2u2R12+u22×p2(β,u2,0)×sgn(sin(β+tan-1u2R1-θ))(R1(u2+h2)R12+u22-T1)}|u2=u0,2dβ,
where
(4)β⊥=(-sinβ,cosβ,0),u0,i=R1(x0·β-hi)R1-x0·β⊥,i=1,2,T1=R1(h2+h1)4R12+(h2-h1)2,kε(r)=12ε∫-∞∞(sgn(r′)+1)w(r-r′ε)dr′,
where ε is a small positive number that is determined by the projection overlapping area, and w(r) is a mollification kernel function as follows:
(5)w(r)={CExp(-11-r2),|r|<1,0,|r|≥1,C=(∫-11Exp(-11-r2)dr)-1.(ii) Ifx0 is on an off-middle plane, for obtaining bθ(x0) we need to define the oblique surface which is through the points S1, S2, and x0. Let O′ denote the intersection point between z-axis and the oblique surface, and we know O′=(0,0,v1) from Figure 2 above. Now, we establish the coordinate system and give the relevant parameters in the oblique surface. Let nO′m denote the coordinate system in the oblique surface, where n, m is the same direction with the vectors u1 and O1′S1, respectively. Let ρ,β′ denote a vector from O′ to x0 and an angle variable in the oblique surface, respectively, and R1′ denote the distance from S1 to O′ in the oblique surface.We know that X-ray source is regarded as a fulcrum in FDK algorithm, and an off-middle plane is approximately obtained by inclining the middle plane. Now, based on the idea of FDK algorithm, we give the steps of deduction ofbθ(x0) in the off-middle plane as follows: (i) writing δbθ(x0)=δbθ,1(x0)+δbθ,2(x0) using the variables in the oblique surface; (ii) finding the relation between δβ and δβ′, where δβ is the rotation angle increment in the middle plane, and δβ′ is the rotation angle increment in the oblique plane; (iii) finding the relation between x0 and ρ; (iv) calculating R1′; (v) obtaining bθ(x0) by accumulating all δbθ(x0) for the angle variable β.Without loss of generality, we give the deduction ofbθ,1(x0) referring to the steps above. From the formula (2), making use of the parameters in the oblique surface nO′m, we obtain. 						(6)δbθ,1(x0)=δβ′R1′2-h12(R1′2-h12-ρ·m)2ddu1×{[kε(R1(u1+h1)R12+u12)-kε(R1(u1+h1)R12+u12-T1)]×(R1′2-h12)-h1u1R1′2-h12+u12×p1(β,u1,v1)×sgn(sin(β+tan-1u1R1-θ))R1(u1+h1)R12+u12}|u1=u0,1,v1=v0,1,
where v0,1=(R1×z0)/(R1-x0×β⊥). Obviously, ρ is x0, when z0=0.We can obtain the relation betweenδβ and δβ′ from Figure 2 as
(7)δβ′=δβR12+h12R12+h12+v12,
and the relation between x0 and ρ as
(8)x0=ρ+v1z.
Since ρ is in the plane nO′m, we can obtain
(9)ρ·(n×m)=0,ρ·m=R1′R12+h12x0·β⊥.
We easily calculate the distance R1′ from S1 to the detectors in the oblique surface,
(10)R1′=|S1O′|=R12+h12+v12.
Now, we substitute formulas (7), (8), (9), and (10) to (6) and obtain bθ,1(x0) by accumulating all δbθ(x0) for the angle variable β as
(11)bθ,1(x0)=∫02π‍U1ddu1×{[kε(R1(u1+h1)R12+u12)-kε(R1(u1+h1)R12+u12-T1)]×(R12+v12)-h1u1R12+v12+u12×p1(β,u1,v1)×sgn(sin(u1R1β+tan-1u1R1-θ))R1(u1+h1)R12+u12}|u1=u0,1,v1=v0,1dβ,
where
(12)U1=(R12+h12)3/2R12+v12M1R12+h12+v12,M1=((R12+v12)(R12+h12)-R12+h12+v12x0·β⊥)2.
Similarly, we can obtain
(13)bθ,2(x0)=∫02πU2ddu2×{(kε(R1(u2+h2)R12+u22-T1))×(R12+v22)-h2u2R12+v22+u22×p2(β,u2,v2)×sgn(sin(β+tan-1u2R1-θ))R1(u2+h2)R12+u22-T1}|u2=u0,2,v2=v0,2dβ,
where
(14)U2=(R12+h22)3/2R12+v22M2R12+h22+v22,M2=((R12+v22)(R12+h22)-R12+h22+v22x0·β⊥)2,
where v0,i=(R1×z0)/(R1-x0×β⊥). When z0=0, the formulas (11) and (13) are (2) and (3), respectively.We need two steps for reconstructingf(x). Firstly, we need to obtain the Hilbert image Hθf(x) of each layer according to the reference [8]. Using the relation of the Hilbert image and the DBP image, we can obtain from the formulas (1), (2), (3), (11), and (13),
(15)Hθf(x)|z=z0=-bθ(x)|z=z02π,
where x=(x,y,z) is an arbitrary point on the inspected object. Secondly, making use of the virtual trajectories and virtual PI-lines in [12, 13], we can obtain the following formula based on the finite inversion formula of Hilbert transform:
(16)f(x)|z=z0=-1(x·θ′⊥-Lt,z)(Ut,z-x·θ′⊥)×(∫Lt,zUt,z(s-Lt,z)(Ut,z-s)×Hθf((x·θ′)θ′+sθ′⊥+z(0,0,1))π(x·θ′⊥-s)ds+Ct,z{∫Lt,zUt,z(s-Lt,z)(Ut,z-s)})|z=z0,
where θ′=(cosθ,sinθ,0), θ′⊥=(-sinθ,cosθ,0), and the constants Lt,z, Ut,z, and Ct,z relate to x·θ′=t and z. We can obtain Ct,z from the integral of f(x) along the virtual PI-lines [14]. However, since the PI-lines are virtual except the middle plane, we cannot obtain an accurate Ct,z. So the proposed reconstruction formula is approximated on the off-middle planes.


## 4. Numerical Experiments

To validate our algorithm, we perform some numerical experiments with simulated data in this section. 3D Shepp-Logan phantom is used to get the simulated data. The parameters of the phantom in 3D Cartesian coordinate systemOxyz are listed in Table 1, where (xc,yc,zc) is the center coordinate of an ellipsoid, and three variables a,b and c denote the length of three half axis of an ellipsoid, and α is the rotating angle of an ellipsoid in the plane xOy. The parameters of two pairs of X-ray sources and panel detectors in the scanning geometry are the same as follows: R1=1426.0mm, R2=1800.0mm, and the length and width of the panel detectors L=W=77.1mm, which is composed of 257 × 257 detector cells.Table 1
Parameters of 3D Shepp-Logan phantom.
N
x
c
y
c
z
c
a
b
c
α (degree)
Density value
1
0
0
0
97.9
73.4
26.8
0
1000.0
2
0
0
0
93.0
70.5
26.2
0
300.0
3
−23.4
0
0
17.0
43.6
20.8
108
-
200
.
0
4
23.4
0
0
11.7
33.0
20.8
72
-
200
.
0
5
0
37.2
0
26.6
22.4
14.9
0
200.0
6
0
10.6
0
4.9
4.9
1.4
0
200.0
7
−6.4
−69.2
0
2.5
4.9
0.6
0
200.0
8
6.4
−69.2
0
2.5
4.9
0.6
90
200.0
9
8.5
0
18.6
5.9
4.3
6.0
0
200.0
10
0
0
18.6
5.9
5.9
6.0
0
200.0According to the parameters above, we can calculate the horizontal diameter of the maximum FOV formed by a set of X-ray source and panel detector is 60.82 mm. From the parameters of 3D Shepp-Logan phantom in Table1, we calculate the longest axis of ellipsoid along x-axis to be 97.9 mm which is greater than 60.82 mm. So we cannot reconstruct its CT image using Cardiac dual-source CT scanning geometry. For acquiring the complete CT data of the 3D Shepp-Logan phantom, we adapt the nonconcyclic dual-source circular cone beam scanning geometry with a larger FOV in Section 2. In this scanning geometry, we set h1=30.42mm and h2=91.26mm which satisfie the conditions above. Each panel detector takes 720 projections within the angle range from 0 to 2π. Two digital radiography (DR) images under the 100th projection angle are shown in Figure 3. In CT image reconstruction, we use θ=0 as the direction of Hilbert transform, and ε=0.24. We reconstruct the CT images from the formulas (15) and (16). Three-image matrix of a horizontal slice is 701×701, as shown in Figure 4. Two-image matrix of a slice perpendicular to x-axis is 701×200, as shown in Figure 5. Two-image matrix of a slice perpendicular to y-axis is 701×200, as shown in Figure 6. For two three-dimensional projection data whose size is 720×257×257, CT image reconstruction takes about 1200 seconds by using CPU and 17.778 seconds by using GPU, where the size of CT image is 701×701×200.Two DR images under the 100th project angle, (a) DR image from the cone beam formed by the X-ray sourceS1 and the panel detectors D1; (b) DR image from the cone beam formed by the X-ray S2 and the panel detectors D2.
(a)
(b)Three reconstruction images of the horizontal slices, (a)z=-16.8mm, (b) z=0, and (c) z=12.8mm.
(a)
(b)
(c)Two reconstruction images of the slices perpendicular tox-axis, (a) x=-20.4mm, (b) x=2.1mm.
(a)
(b)Two reconstruction images of the slices perpendicular toy-axis, (a) y=-6.0mm, (b) y=0.9mm.
(a)
(b)


## 5.  Conclusion

In this work we give nonconcyclic dual-source circular cone beam scan mode with a larger FOV, and deduce the reconstruction formula. Using two pairs of projections obtained from this nonconcyclic dual-source scanning geometry, the FOV is enlarged effectively in the same equipment condition without data rebinning. It is because of the fact that for the real CT system, the flux output from X-ray source is not isotropic, and then the data acquired in this scanning mode using small panel detectors are relative to more uniform intensity of the flux output than a large one. The experiment confirmed that our reconstruction method is effective, when a cone angle is small.



---

*Source: 526247-2013-03-19.xml*