Exact fan-beam image reconstruction algorithm for ... - Semantic Scholar
INSTITUTE OF PHYSICS PUBLISHING Phys. Med. Biol. 50 (2005) 1805–1820
PHYSICS IN MEDICINE AND BIOLOGY
doi:10.1088/0031-9155/50/8/012
Exact fan-beam image reconstruction algorithm for truncated projection data acquired from an asymmetric half-size detector Shuai Leng1, Tingliang Zhuang1, Brian E Nett1 and Guang-Hong Chen1,2 1 2
Department of Medical Physics, University of Wisconsin–Madison, WI 53704, USA Department of Radiology, University of Wisconsin–Madison, WI 53792, USA
E-mail: [email protected]
Received 19 November 2004, in final form 7 February 2005 Published 6 April 2005 Online at stacks.iop.org/PMB/50/1805 Abstract In this paper, we present a new algorithm designed for a specific data truncation problem in fan-beam CT. We consider a scanning configuration in which the fan-beam projection data are acquired from an asymmetrically positioned half-sized detector. Namely, the asymmetric detector only covers one half of the scanning field of view. Thus, the acquired fan-beam projection data are truncated at every view angle. If an explicit data rebinning process is not invoked, this data acquisition configuration will reek havoc on many known fan-beam image reconstruction schemes including the standard filtered backprojection (FBP) algorithm and the super-short-scan FBP reconstruction algorithms. However, we demonstrate that a recently developed fan-beam image reconstruction algorithm which reconstructs an image via filtering a backprojection image of differentiated projection data (FBPD) survives the above fan-beam data truncation problem. Namely, we may exactly reconstruct the whole image object using the truncated data acquired in a full scan mode (2π angular range). We may also exactly reconstruct a small region of interest (ROI) using the truncated projection data acquired in a short-scan mode (less than 2π angular range). The most important characteristic of the proposed reconstruction scheme is that an explicit data rebinning process is not introduced. Numerical simulations were conducted to validate the new reconstruction algorithm. (Some figures in this article are in colour only in the electronic version)
1. Introduction Recently, developments of novel and mathematically exact fan-beam image reconstruction methods have attracted much attention and many discoveries have been made. These major 0031-9155/05/081805+16$30.00 © 2005 IOP Publishing Ltd Printed in the UK
1805
1806
S Leng et al
breakthroughs address two important issues in fan-beam CT, namely, the possibility of a mathematically exact reconstruction of a region of interest (ROI) in the super-short scan mode, and reconstruction from partially truncated projection data. The first issue is the discovery of the super-short scan mode and corresponding fan-beam image reconstruction algorithms (Noo et al 2002, Kudo et al 2002, Chen 2003, Pan et al 2004, Zhuang et al 2004). These developments are based on either the Radon transform and a novel data rebinning scheme (Noo et al 2002, Kudo et al 2002), a Tuy-like intermediate function (Chen 2003, Zhuang et al 2004), or a dimension reduction scheme from cone-beam conventional helical CT to fan-beam CT (Pan et al 2004). These new algorithms may also be divided into two classes. The first class (Noo et al 2002, Kudo et al 2002, Chen 2003) may be viewed as a fan-beam counterpart of the Katsevich algorithm for the conventional helical CT (Katsevich 2002a, 2002b, 2004). The second class (Pan et al 2004, Zhuang et al 2004, Noo et al 2004) may be viewed as a fan-beam counterpart of the Zou–Pan algorithms for the conventional helical CT (Zou and Pan 2004a, 2004b). The second aspect is to develop new methods of image reconstruction from partially truncated projection data (Clackdoyle et al 2004, Noo et al 2004). When the novel concept of a virtual scanning trajectory is introduced (Clackdoyle et al 2004), it has been demonstrated that the first class of super-short scan reconstruction algorithms (Noo et al 2002, Kudo et al 2002, Chen 2003) may be used to reconstruct images from a partially truncated projection data set. In this method, it seems that the data rebinning process is indispensable. This rebinning algorithm addresses the problem in which the projection data have been truncated only for some view angles. However, it has also been demonstrated that the partial data truncation problem may be solved by using the second class of the reconstruction algorithms (Noo et al 2004). In this case, it is not necessary to explicitly introduce the data rebinning process. A feature of this truncated scanning configuration is that projection data may or may not be truncated for each view angle. The necessary condition for exact reconstruction of an ROI is that the projection data corresponding to the rays passing through the ROI are untruncated in all views. In addition, a large class of new image reconstruction formulae for fan-beam projection data have also been discovered, and have been proven to be equivalent to Radon’s original formula (Clackdoyle and Noo 2004). In this paper, we consider a specific fan-beam scanning configuration with data truncation at every view angle: fan-beam projection data are acquired using an asymmetric detector that covers only one half of the image object (figure 1). We refer to this scenario as a constant view truncation (CVT) problem. This is a special case of the more general partial view data truncation problems that frequently occur in nuclear medicine (Chang et al 1995, Loncaric et al 1995, La Riviere et al 2001), radiation therapy (Hooper and Fallone 2002), and diagnostic x-ray CT (Hsieh et al 2004), where the detector is large enough to cover more than half of the image object but not large enough to cover the whole image object at some view angles. If the data are acquired from a complete circle, it is known that the data set is sufficient to reconstruct the whole image object (Chang et al 1995). However, the standard fan-beam filtered backprojection (FBP) image reconstruction algorithm does not provide a mathematically exact image reconstruction without an explicit data rebinning process from truncated fan-beam projections to non-truncated parallel-beam projections or from truncated fan-beam projections to non-truncated fan-beam projections. This new algorithm provides a means for exact reconstruction of the entire image object, without invoking a rebinning operation, from truncated data. In this paper, we will demonstrate that the newly developed fan-beam image reconstruction algorithm (Zhuang et al 2004) via filtering the backprojection image of differentiated projection data (FBPD) provides a mathematically exact reconstruction algorithm for the aforementioned CVT problem. Using the FBPD algorithm, data rebinning is no longer
Exact fan-beam image reconstruction algorithm
Focal Spot Trajectory
1807
right side
left side
Region of Image Support
(Co
nst Half ant Det Vie etec w T tor run ca
tion
) Full Detector
(No Truncation)
Figure 1. The specific data truncation problem which will be addressed in this paper is that in which projection data are collected from a detector of only half the size required for a non-truncated data. The detector is positioned asymmetrically for every given view angle. Note that the definitions of the left-hand side of the detector and the right-hand side of the detector are taken from the source position looking towards the detector.
needed. On the basis of an appropriate choice of filtering lines, a local ROI may be exactly reconstructed using the truncated projection data acquired from a scanning path which is less than 2π in angular coverage. 2. Brief review of the fan-beam FBPD reconstruction method: without data truncation For simplicity, we consider image reconstruction from fan-beam projections on an arc source trajectory with radius r (figure 2). The source trajectory is parametrized as y(t) = r(cos t, sin t),
t ∈ [ti , tf ],
(2.1)
where t parametrizes the view angle and is measured counter-clockwise with respect to the x-axis, ti and tf are the view angles corresponding to the starting and ending points of the scanning path. The newly developed FBPD image reconstruction scheme (Zhuang et al 2004) can be briefly summarized as below. For a given image point x, suppose a straight line, l, x )) that passes through the image point x intersects the source trajectory at two points y(t a ( x )). According to the Tuy data sufficiency condition (Chen 2003, Tuy 1983), the and y(t b ( x ), tb ( x )] are sufficient projection data from the scanning path within the angular range [ta ( for the reconstruction of the image point x. Throughout this section, we assume that only x ), tb ( x )] be used. This amounts to a special the projection data within the angular range [ta (
1808
S Leng et al
Figure 2. Definition of each of the vectors and angles used in deriving the FBPD reconstruction formula.
‘one or zero’ weighting scheme among a large class of other potential weighting schemes (Zhuang et al 2004). On the basis of this specific choice of the weighting function, the filtering must be performed along the straight line l in the space of the backprojection image x ), tb ( x )], one may construct a of the differentiated data. Using the projection data from [ta ( backprojection image Q( x , x ), tb (x ) ∂ ∂ 1 − [gm (γ = γx , t) − gm (γ = γx + π, t)], (2.2) Q( x , x ) = dt | x − y (t)| ∂t ∂γ ta ( x) where gm (γ , t) is the measured fan-beam projection data at view angle t and fan angle γ , where γ is measured from the iso-ray using the counter-clockwise convention. The values of γx and βˆ x are determined by (figure 2) γx = φx − t − π,
βˆ x = (cos φx , sin φx ),
(2.3)
where x is an arbitrary point in the backprojection image space. After the formation of a backprojection image associated with the image point x, the final image reconstruction is obtained through a Hilbert transform along the filtering line l. Namely, Q( x , x )|l 1 f ( x ) = − 2 dxl , (2.4) 2π l xl − xl where xl and xl are the projection of vectors x and x onto the filtering line l. For clarity, one may choose a coordinate system such that the x-axis is parallel to the filtering line l. In this case the above filtering process may be explicitly written as +∞ Q(x, y; x , y) 1 dx . (2.5) f (x, y) = − 2 2π −∞ x − x Note that the range of integration extends to infinity. Thus, we need to calculate the backprojection image values over an infinite range. This is mathematically possible as demonstrated in the above formulae (2.2) and (2.3). However, it is not computationally tolerable. Fortunately, for a compactly supported image function f ( x ) the image values of f ( x ) on the nonzero portion of the filtering line may be exactly calculated using a finite range Hilbert transform (Tricomi 1957, Milkhlin 1957). If we assume that the image to be reconstructed is nonzero only within a circle of radius R0 , and the filtering operation is
Exact fan-beam image reconstruction algorithm
1809
performed along a horizontal line (l) as indicated by equation (2.5), then the following formula is suggested to reconstruct image points ( x ): √ +R 1 R 2 − x 2 1 Q( x , x ) + C , dx (2.6) f (x, y0 ) = − 2 √ 2π x − x R 2 − x 2 −R where x = (x, y0 ), x = (x , y0 ), R R0 and C is a constant which is filtering line dependent. x , x ) are Using this formula, only the data in the range x ∈ [−R, +R] of the function Q( needed in the image reconstruction. This finite Hilbert transform reduces the computational load significantly since we now only need to calculate the values of the backprojection image within a finite range. To determine the constant C in equation (2.6), one can use the fact that the image values should be zero when the image point is outside the region of support (Zhuang et al 2004, ¯ y0 ) is outside the image support, then the Noo et al 2004). Therefore, if a point xc = (x, constant C can be determined as √ +R R 2 − x 2 Q( x c , x ). dx (2.7) C=− −x ¯ x −R In practice, several points along the filtering lines, which were also outside the region of compact support, were chosen to calculate the constant C. The mean value of several calculated C values was used in (2.6) to determine the image values for points along the same filtering line. An alternative means to determine the constant C is to integrate the both sides of equation (2.6) along the filtering line l. Namely, integration over the variable x is calculated. The integral on the left-hand side gives nothing but the projection value along the filtering line and it is denoted as gm (γ0 , t0 ). The integral over the variable x in the first term of the right-hand side of equation (2.6) vanishes. The integral of the second term is given by −C/2π . Thus, the constant C is determined by the following equation (Tricomi 1957): C = −2πgm (γ0 , t0 ).
(2.8)
The advantage of this method to determine the constant C is that there is no need to calculate the image values for the points outside the image support. Therefore, we only need the backprojection image values in the range [−R0 , +R0 ]. In summary, we emphasize the following important feature of the FBPD method given by equations (2.2) and (2.6): if the value of the backprojection image Q( x , x ) is correct for all points along a filtering line then any point on that filtering line may be accurately reconstructed. Therefore, it is crucial to obtain correct values of the backprojection image Q( x , x ) for all the points along the filtering line l. The correct values of the backprojection image Q( x , x ) x ), tb ( x )] provides a projection value for any are achieved when: every view angle t ∈ [ta ( given point x on a chosen filtering segment. In the case of non-truncated data as discussed in Clackdoyle and Noo (2004), in order to obtain the correct value of the backprojection image Q( x , x ), there is a degree of freedom in the selection of the scanning path of the backprojection region for a given filtering line. Namely, one may backproject the differentiated data from either of two given arcs (figure 3). One may also selectively combine the projection values of the conjugate rays from two arcs in figure 3 to form the correct backprojection image. The use of compensation from conjugate rays will be elaborated upon below. 3. Fan-beam FBPD image reconstruction: constant view truncation (CVT) In the standard fan-beam filtered backprojection (FBP) reconstruction (Kak and Slaney 1987) algorithm and the super-short scan image reconstruction algorithms developed in Noo et al
1810
S Leng et al
ROI
ROI
Filtering Line (l)
(a)
Filtering Line (l)
(b)
Figure 3. Options of angular range for a backprojection image from non-truncated projection data. In order to achieve the correct backprojection image for all points on the filtering line given above, one may select the data from either the source trajectory highlighted in (a) or the trajectory highlighted in (b).
(2002), Kudo et al (2002) and Chen (2003), the filtering step is conducted over all of the projection data for each given view angle. The filtered data are then backprojected in order to form the reconstructed image. Note that the filtering kernels, either a ramp filter or a Hilbert filter, in the aforementioned FBP-type image reconstruction algorithms are non-local. Therefore, for the truncated projection data, an explicit rebinning process is necessary for the use of the previously mentioned FBP image reconstruction algorithms. Originally, the FBPD image reconstruction algorithm (Zhuang et al 2004) was derived under the assumption that the projection data are not truncated for any view angle. In this section, we will elaborate on this point carefully to demonstrate that the final FBPD image reconstruction algorithm given by equations (2.2) and (2.6) survives the CVT problem. In order to address the CVT problem using the FBPD algorithm, the following observations are crucial: (1) for a given point x in the backprojection image space, only one projection for each view angle is required to construct the backprojection image. Therefore, for a given x ), tb ( x )] and a given point x in the backprojection image space, one view angle t ∈ [ta ( does not have to know all of the projection values within the fan-shape region that covers the whole image object in order to properly reconstruct the image. (2) For a given point x in the backprojection image space, only one ray is required in (1). This ray may be from a view x ), tb ( x )] or it may be a conjugate ray from a view angle out of the angular range angle t ∈ [ta ( x ), tb ( x )]. The second observation implies that one may combine the projection data from [ta ( different view angles in order to form a correct backprojection image. We analyse below how to apply the FBPD reconstruction algorithm to the CVT problem (figure 1), where an asymmetrically placed half-size detector is utilized to acquire fan-beam projection data. In order to facilitate the analysis, we assume that the image function has a definite boundary (figure 4). Outside the boundary , the image function f ( x ) ≡ 0. For a given image point x, the filtering line l intersects the boundary curve at a minimum of two points A and B (figure 4). The detector is large enough to cover one half of the image object. Namely, if the iso-ray divides the image object into two parts at any view angle, the detector is always sufficiently large to guarantee that there are no data truncation for one of the two parts. Each of the schematic diagrams in the present paper assumes the convention that the data provided are on the right-hand side of the iso-ray given that one is oriented from the source towards the detector (figure 1). A third-generation scanning configuration and an equal-angular detector
Exact fan-beam image reconstruction algorithm
l A
x
B
1811
A
B
l
x
Γ
Γ
(a)
(b)
Figure 4. Examples of possible filtering lines l, and the definition of points A and B which are the outer most points intersection points between the filtering line and the curve . Source Trajectory Illuminating Sample Image Point
F2
F1
F3
x x
(a)
x
x
(b)
(c)
(d)
Figure 5. Demonstration that in the CVT case any one given image point is only illuminated by rays originating from an arc of arc-length π . Three sample focal spots F1 , F2 and F3 are given to illustrate this point in (a), (b) and (c); the total angular range is shown in (d).
are also assumed. The analysis may easily be adapted to the case of a collinear detector or the case of a left-hand side convention of object illumination. Using the above scanning configuration, for a given image point x, it is easy to see that only one half of the circular scanning path will illuminate the image point x (figure 5). Geometrically, if one draws a straight line that passes through the points x and the iso-centre, the straight line bisects the circular scanning path. According to the convention of the scanning configuration, one may easily determine which half of the scanning path will illuminate the point x. We are now in the position to discuss the impact of the selection of the filtering line on the reconstruction algorithm. The choice of filtering lines in the reconstruction algorithm will determine the angular range of data required to achieve an exact reconstruction for all of the points on that given line. In order to provide a general foundation for the selection of filtering lines two distinct classes of filtering lines will be discussed below. The first class is all the filtering lines which pass through the origin, and the second class is the filtering lines which do not pass through the origin. The angular range of data required for each of these scanning paths will also be addressed. Following this analysis will be a discussion of how one may reconstruct an ROI from a group of filtering lines. In the first class of filtering lines, one selects a filtering line as the line that connects an image point x and the iso-centre (figure 4(a)). Note that one half of the scanning path will illuminate boundary point A. In order to illuminate the other boundary point B, data are required from the other half of the scanning path. Therefore, only a full scan can generate
1812
S Leng et al Required Source Trajectory
ROI
Required Source Trajectory ROI
Correct Backprojection Data
ROI
Filtering Line (l)
Correct Backprojection Data
Correct Backprojection Data
Required Source Trajectory
(a)
(b)
(c)
Figure 6. In comparison with the case of untruncated data in figure 3 in order to provide the correct backprojection data for a central filtering line a full scan is needed for the case of constant view truncation. In this case it is evident that the union of the backprojection data required in (a) and (b) form the complete 2π trajectory (c).
sufficient projection data for all the points along the radial line (figure 6). This point may be further clarified by analysing the illumination pattern given in figures 5(a)–(c). When the focal spot moves from figure 5(a) to figure 5(c), only the left-hand side of the horizontal line is illuminated. In other words, one may obtain the correct values of the backprojection image for all the points on the left portion of the horizontal filtering line (figures 5(d) and 6(a)). In order to obtain the correct backprojection image values required for the right portion of the horizontal filtering line, one must scan through the other half of the circular scanning path (figure 6(b)). Therefore, a full 2π angular range is required to provide sufficient data to generate the correct backprojection images for all of the points lying along this filtering line (figure 6(c)). This situation is significantly different from the non-truncated case illustrated in figure 3. Due to the symmetry breaking when moving from the non-truncated data reconstruction problem to the half-size detector CVT reconstruction problem, the degree of freedom in the backprojection step to choose either arc (figure 3(a) or figure 3(b)) has been lost. In the CVT reconstruction case in order to generate correct values for the backprojection image at each point along this radial filtering line data are required from the full angular range (figure 6(c)). The second class of filtering lines is composed of those filtering lines that do not pass through the iso-centre. One such example is illustrated in figure 7. In this case, two half circle illumination regions corresponding to the boundary points A and B overlap with one another. The union of these illumination regions is less than a full circle. In the radial filtering case discussed above the union of these two illumination regions consisted of the full 2π trajectory (figure 6). The starting and ending view angles will be denoted as ti and tf . It is important to observe that the view angles in the angular range [ta , ti ], which are formally required by the FBPD reconstruction algorithm, are missing. Thus, there are no rays emanating from this angular range. The data in the angular range [ta , ti ], which are required in order to generate the correct values of the backprojection image associated with point A are not present in the region defined by the union of the illumination regions for points A and B. However, one can use the conjugate rays emanating from the view angles in the angular range [ti , tf ] to compensate for the missing projections in the angular range [ta , ti ] (figure 8). Therefore, the angular range [ti , tf ] is sufficient for the reconstruction of the filtering segment AB . To briefly summarize the analysis given above, we conclude that the selection of the filtering lines is an important issue and has implications on the required angular range of projection data. Now that the effect of the orientation of the filtering lines has been probed we may address the issue of reconstructing an ROI from a combination of filtering lines. In the case of radial
Exact fan-beam image reconstruction algorithm
1813
l
A
B
x
Γ
Figure 7. In order to reconstruct x and all other points along the filtering line l the required source trajectory is given by the union of the all source positions that illuminate image point A and all the source positions that illuminate image point B. Conjugate Ray Conjugate Source Position
l
Missing Source Position
ti ta
Conjugate Rays
B
A
Γ
l
ti ta
B
A
Γ
tf
tf (a)
(b)
Figure 8. The data in the angular range [ta , ti ] are missing. However, since for each view angle the backprojection image only has a contribution from a single projection, the missing projections may be compensated for by substituting their equivalent conjugate projections. The conjugate rays are emitted from different view angles (a). In each case the conjugate source point is within the union of the the illumination regions of points A and B (b). Remember that the right-hand convention is still used here and thus it is obvious that the conjugate ray will not be truncated.
filtering (figure 4(a)), there is no small ROI reconstruction. Since the union of the regions of illumination for any given radial filtering lines composes the entire 2π circle (figure 6), there is no possibility of a ‘short scan’. However, in the case of non-radial filtering lines (figure 4(b)) the possibility of a ‘short scan’ exits. As shown in figure 9 the angular range required for the outermost filtering line of the scanning path is sufficient to provide a complete backprojection data set for an ROI bounded by the filtering line and the boundary curve . Thus, for non-radial filtering lines one may accurately reconstruct a local ROI using a ‘short scan’ mode. Namely, the angular range of
1814
S Leng et al Required Source Trajectory
Reconstructed ROI B
l
B
B
l
A
B
l l
A
A
Γ
(a)
Γ
(b)
Γ
(c)
A
Γ
(d)
Figure 9. Demonstration that the angular range required to generate the backprojection data for the outermost filtering line (a) is sufficient to generate the backprojection data for an ROI. This illustrative case shows that an ROI may be reconstructed using filtering lines parallel to the outermost filtering line (b) and (c). Thus, the FBPD provides the opportunity for a ‘short scan’ from CVT data (d).
the scanning path is less than 2π . Note that the meaning of our short scan is different from the conventional case with scanning angular range π + γm (γm is the fan-angle of the whole image object) (Kak and Slaney 1987). The selection of the optimal filtering lines is dependent upon the given scanning configuration and the desired ROI, and thus the topic of the selection of the optimal filtering lines will not be addressed here. It is important to emphasize that the short-scan mode is a significant result of the FBPD based image reconstruction algorithm. In the conventional strategy of rebinning the truncated fan-beam projections to non-truncated parallel beam projections, no short-scan mode is possible. Namely, even for the case of small ROI reconstruction, data from the entire scanning path are required. 4. Numerical simulations In this section, we present numerical simulations in order to illustrate the validity of the proposed fan-beam FBPD image reconstruction from the constantly truncated projection data. A standard Shepp–Logan phantom (Kak and Slaney 1987) was used to generate the projection data from an arc source trajectory. The Shepp–Logan dimensions were scaled by a constant factor of 5.8. The total fan angle was chosen to be γm = π4 . The sampling rate of the equiγm . The sampling rate of the view angles was chosen angular detector was chosen as γ = 723 π to be t = 1031 . Each of the images displayed in this section has a compressed window of [0.98 1.06] in order to highlight the low contrast objects in the phantom. In computing the backprojection image of differentiated projection data, the derivatives were calculated by a three-point formula. A Hann window was used to suppress the highfrequency fluctuations in the reconstructed images. The constant C in equation (2.6) was determined by using equation (2.8). A bilinear interpolation scheme was used to find the projection values for a given filtering line. Figure 10 shows a simulation for a complete scan, where horizontal parallel filtering lines were chosen. Figure 10(a) demonstrates the scanning path and orientation of the filtering lines. Figure 10(b) shows the reconstructed image. Figures 10(c) and 10(d) present the central horizontal line and the central vertical line of the image in figure 10(b). Theoretical values are presented for comparison in each case. The simulation results show that the whole object can be reconstructed correctly when a full scan is used with the asymmetric half detector.
Exact fan-beam image reconstruction algorithm
A
1815
B
Full FOV Supported (a) 1.08
(b)
Theoretical Values Reconstructed Values
Theoretical Values Reconstructed Values
1.06
1.06 1.05 1.04 1.04 1.02 1.03 1
1.02
0.98
1.01
0.96
1 –4
–3
–2
–1
0
(c)
1
2
3
4
–6
–4
–2
0
2
4
6
(d)
Figure 10. Simulation results of a full scan (2π ) source trajectory. A schematic figure illustrating the scanning path used and the orientation of filtering lines (a). Parallel horizontal filtering lines were used here. The reconstructed image (b) along with density plots of the central horizontal and vertical lines (c) and (d) demonstrate good agreement with theoretically expected values.
Figures 11 and 12 are simulation results of reconstructing a specific ROI using a short scan. The scanning path was parametrized as (2π − α), where α = π/4 for each of the ROI simulations given below. Figures 11(a) and 12(a) demonstrate the scanning path and direction of the filtering lines for ROIs in the upper and lower portion of the image respectively. Figures 11(b) and 12(b) demonstrate the reconstructed images. Figures 11(c) and 12(c) present density plots from horizontal lines inside of ROI, and 11(d) and 12(d) present horizontal lines outside of ROI. The density values within the ROI have been accurately reconstructed. Quantitatively exact CT number cannot be obtained outside the ROI, but some of the image structure remains visible in regions near the ROI. Figures 11(e) and 12(e) plot the density values along central vertical lines of the reconstructed images. The region of support for exact reconstruction is defined by (y > 1.6) and (y < −1.6) for figures 11 and 12 respectively. In each case the theoretical phantom values are presented in each plot for comparison. Figure 13 is the result of a short scan where the ROI is the left most portion of the phantom. The main difference between this reconstruction and that presented in
1816
S Leng et al
Supported ROI B
A
γm
(a)
(b) 1.08
Theoretical Values Reconstructed Values
1.08
Theoretical Values Reconstructed Values
1.06 1.06 1.04 1.04 1.02
1.02
1
1
0.98
0.98
0.96
0.96 –4
–3
–2
–1
0
1
2
3
4
–4
(c)
–3
–2
–1
0
1
2
3
4
(d)
Theoretical Values Reconstructed Values
1.06
1.05
1.04
1.03
1.02
1.01
1
–6
–4
–2
0
2
4
6
(e)
Figure 11. Reconstruction of an ROI from the top part of the object using a short-scan source trajectory. A schematic figure of this short-scan trajectory and orientation of the outermost filtering line is given in (a). Parallel horizontal filtering lines were used to reconstruct the ROI. The dashed arc presents the angular range of source trajectory. Reconstructed image where the ROI is given by the region above the solid line (b). Density plot of two horizontal lines (c) and (d) at the positions indicated by the dashed lines in (b), where (c) is within the ROI and (d) is outside the ROI. The density plot of the central vertical line is given in (e), where the image values for the exact ROI correspond to y > 1.6.
Exact fan-beam image reconstruction algorithm
1817
γm A
B
Supported ROI (a)
(b)
Theoretical Values Reconstructed Values
1.08
1.08
1.06
1.06
1.04
1.04
1.02
1.02
1
1
0.98
0.98
Theoretical Values Reconstructed Values
0.96
0.96 –4
–3
–2
–1
0
1
2
3
4
(c) 1.06
–4
–3
–2
–1
0
1
2
3
4
(d)
Theoretical Values Reconstructed Values
1.05
1.04
1.03
1.02
1.01
1
–6
–4
–2
0
2
4
6
(e)
Figure 12. Reconstruction of an ROI from the bottom portion of the phantom using a short-scan source trajectory. A schematic figure of this short-scan trajectory and orientation of the outermost filtering line is given in (a). Parallel horizontal filtering lines were used to reconstruct the ROI. The dashed arc presents the angular range of source trajectory. Reconstructed image where the ROI is given by the region below the solid line (b). Density plot of two horizontal lines (c) and (d) at the positions indicated by the dashed lines in (b), where (c) is within the ROI and (d) is outside the ROI. The density plot of the central vertical line is given in (e), where the image values for the exact ROI correspond to y < −1.6.
1818
S Leng et al
A
Supported ROI
B
γm
(a) 1.08
(b)
Theoretical Values Reconstructed Values
Theoretical Values Reconstructed Values
1.08
1.06 1.06 1.04 1.04 1.02 1.02 1 1 0.98 0.98 0.96
0.96 –4
–3
–2
–1
0
1
2
3
4
(c)
1.08
–5
–4
–3
–2
–1
0
1
2
3
4
5
(d)
Theoretical Values Reconstructed Values
1.06
1.04
1.02
1
0.98
0.96 –4
–3
–2
–1
0
1
2
3
4
(e)
Figure 13. Reconstruction of an ROI from the left portion of the object using a short-scan source trajectory. A schematic figure of this short-scan trajectory and orientation of the outermost filtering line is given in (a). Parallel vertical filtering lines were used to reconstruct the ROI. The dashed arc presents the angular range of source trajectory. Reconstructed image where the ROI is given by the region to the left of the solid line (b). Density plot of two vertical lines (c) and (d) at the positions indicated by the dashed lines in (b), where (c) is within the ROI and (d) is outside the ROI. The density plot of the central horizontal line is given in (e), where the image values for the exact ROI correspond to x < −1.6.
Exact fan-beam image reconstruction algorithm
1819
figures 11 and 12 is that in this case vertical filtering lines were chosen. Figure 13(a) demonstrates the scanning path and orientation of the filtering lines. Figure 13(b) demonstrates the reconstructed image. Figures 13(c) and (d) present vertical lines which lie both inside and outside of the ROI. Figure 13(e) presents the density values of the central horizontal line in the reconstructed image. Note that in figures 11(e), 12(e) and 13(e) as expected the reconstructed values match the theoretical values within each ROI, and as the distance increases from the ROI the reconstructed values become increasingly less accurate. In summary fan-beam image reconstruction simulations were performed using the FBPD reconstruction algorithm on CVT data for scans of both the entire object as well as specific ROI reconstructions. In the case of ROI reconstructions a ‘short scan’ was used in which the angular range was less than 2π . The simulations support the theoretical assertions made above, and demonstrate that the FBPD algorithm provides accurate reconstruction of the entire object (full-scan case) or of a given region of interest (short-scan case). 5. Conclusions and discussion In this work, we present a new image reconstruction algorithm for a specific data truncation problem. The scanning configuration is that of a third generation geometry CT scanner in which only a detector of half of the size required for non-truncated projections is placed asymmetrically about the iso-ray. In this geometry, the data are truncated at every view angle, and thus this problem has been referred to as constant view truncation. If the scanning path is a complete circle, it was previously demonstrated that the data are sufficient to reconstruct the whole image object (Chang et al 1995) after the fan-beam projection data are rebinned into parallel-beam projection data. In this rebinning scheme, even in the case of ROI reconstruction a scanning path shorter than 2π is not sufficient. In this work, we demonstrate that an ROI image may be accurately reconstructed using constantly truncated projection data collected from an angular range of less than the full 2π scan. The connection between the acquired data and the generated backprojection image, in the recently developed FBPD reconstruction (Zhuang et al 2004), has been exploited to yield the feasibility of a ‘short scan’ in the constant view truncation (CVT) reconstruction problem. Numerical simulations were conducted to validate the proposed image reconstruction scheme for both an ROI reconstruction and the reconstruction of a whole object. Acknowledgments This work was partially supported by National Institute of Health grants 1R21 EB001683-01, 1R21 CA109992-01, 5T32CA009206-24, a start-up grant from the University of Wisconsin– Madison, and a grant from GE Healthcare Technologies. References Chang W, Loncaric S, Huang G and Sanpitak P 1995 Asymmetric fan transmission CT on SPECT system Phys. Med. Biol. 40 913–28 Chen G-H 2003 A new framework of image reconstruction from fan beam projections Med. Phys. 30 1151–61 Clackdoyle R and Noo F 2004 A large class of inversion formulae for the 2D Radon transform of functions of compact support Inverse Problems 20 1281–91 Clackdoyle R, Noo F, Guo J and Roberts J 2004 A quantitative reconstruction from truncated projections in classical tomography IEEE Trans. Nucl. Sci. 51 2570–8 Hooper H R and Fallone B G 2002 Sinogram merging to compensate for truncation of projection data in tomotherapy imaing Med. Phys. 29 2548–51
1820
S Leng et al
Hsieh J, Chao E, Thibault J, Grekowicz B, Horst A, McOlash S and Myers T J 2004 A novel reconstruction algorithm to extend the CT scan field-of-view Med. Phys. 31 2385–91 Kak A C and Slaney M 1987 Principles of Computerized Tomographic Imaging (New York: IEEE Press) Katsevich A 2002a Theoretically exact filtered backprojection-type inversion algorithms for spiral CT SIAM J. Appl. Math. 62 2012–26 Katsevich A 2002b Analysis of an exact inversion algorithm for spiral cone-beam CT Phys. Med. Biol. 47 2583–97 Katsevich A 2004 An improved exact filtered backprojection algorithm for spiral computed tomography Adv. Appl. Math. J. Math. Math. Sci. 32 681–97 Kudo H, Noo F, Defrise M and Clackdole R 2002 New super-short scan algorithms for fan-beam and cone-beam reconstruction Conf. record of the 2002 IEEE Nuclear Science Symposium and Medical Imaging Conference (Norfolk, VA) (IEEE Service Centre vol 2) pp 902–6 La Riviere P J, Pan X, Gilland D, Kao C-M, Chang W and Jaszczak R 2001 Transmission image reconstruction and redundant information in SPECT with asymmetric fanbeam collimation IEEE Trans. Nucl. Sci. 48 1357–63 Loncaric S, Chang W and Huang G 1995 A processing technique for the truncated projections of asymmetric-fan-beam transmission imaging IEEE Trans. Nuc. Sci. 42 2292–7 Milkhlin S G 1957 Integral Equations (New York: Pergamon) pp 126–31 Noo F, Defrise M, Clackdoyle R and Kudo H 2002 Image reconstruction from fan-beam projections on less than a short-scan Phys. Med. Biol. 47 2525–46 Noo F, Clackdoyle R and Pack J D 2004 A two-step Hilbert transform method for 2D image recontruction Phys. Med. Biol. 49 3903–23 Pan X, Xia D, Zou Y and Yu L 2004 A unified analysis of FBP-based algorithms in helical cone-beam and circular cone- and fan-beam scans Phys. Med. Biol. 49 4349–69 Tricomi F G 1957 Integral Equations (New York: Interscience) pp 173–85 Tuy H K 1983 An inversion formula for cone-beam reconstruction SIAM J. Appl. Math. 43 546–52 Zhuang T, Leng S, Nett B E and Chen G-H 2004 Fan-beam and cone-beam image reconstruction via filtering the backprojection image of differentiated projection data Phys. Med. Biol. 49 5489–503 Zou Y and Pan X 2004a Exact image reconstruction on PI-lines from minimum data in helical cone-beam CT Phys. Med. Biol. 49 941–59 Zou Y and Pan X 2004b Image reconstruction on PI-lines by use of filtered backprojection in helical cone-beam CT Phys. Med. Biol. 49 2717–31
PHYSICS IN MEDICINE AND BIOLOGY
doi:10.1088/0031-9155/50/8/012
Exact fan-beam image reconstruction algorithm for truncated projection data acquired from an asymmetric half-size detector Shuai Leng1, Tingliang Zhuang1, Brian E Nett1 and Guang-Hong Chen1,2 1 2
Department of Medical Physics, University of Wisconsin–Madison, WI 53704, USA Department of Radiology, University of Wisconsin–Madison, WI 53792, USA
E-mail: [email protected]
Received 19 November 2004, in final form 7 February 2005 Published 6 April 2005 Online at stacks.iop.org/PMB/50/1805 Abstract In this paper, we present a new algorithm designed for a specific data truncation problem in fan-beam CT. We consider a scanning configuration in which the fan-beam projection data are acquired from an asymmetrically positioned half-sized detector. Namely, the asymmetric detector only covers one half of the scanning field of view. Thus, the acquired fan-beam projection data are truncated at every view angle. If an explicit data rebinning process is not invoked, this data acquisition configuration will reek havoc on many known fan-beam image reconstruction schemes including the standard filtered backprojection (FBP) algorithm and the super-short-scan FBP reconstruction algorithms. However, we demonstrate that a recently developed fan-beam image reconstruction algorithm which reconstructs an image via filtering a backprojection image of differentiated projection data (FBPD) survives the above fan-beam data truncation problem. Namely, we may exactly reconstruct the whole image object using the truncated data acquired in a full scan mode (2π angular range). We may also exactly reconstruct a small region of interest (ROI) using the truncated projection data acquired in a short-scan mode (less than 2π angular range). The most important characteristic of the proposed reconstruction scheme is that an explicit data rebinning process is not introduced. Numerical simulations were conducted to validate the new reconstruction algorithm. (Some figures in this article are in colour only in the electronic version)
1. Introduction Recently, developments of novel and mathematically exact fan-beam image reconstruction methods have attracted much attention and many discoveries have been made. These major 0031-9155/05/081805+16$30.00 © 2005 IOP Publishing Ltd Printed in the UK
1805
1806
S Leng et al
breakthroughs address two important issues in fan-beam CT, namely, the possibility of a mathematically exact reconstruction of a region of interest (ROI) in the super-short scan mode, and reconstruction from partially truncated projection data. The first issue is the discovery of the super-short scan mode and corresponding fan-beam image reconstruction algorithms (Noo et al 2002, Kudo et al 2002, Chen 2003, Pan et al 2004, Zhuang et al 2004). These developments are based on either the Radon transform and a novel data rebinning scheme (Noo et al 2002, Kudo et al 2002), a Tuy-like intermediate function (Chen 2003, Zhuang et al 2004), or a dimension reduction scheme from cone-beam conventional helical CT to fan-beam CT (Pan et al 2004). These new algorithms may also be divided into two classes. The first class (Noo et al 2002, Kudo et al 2002, Chen 2003) may be viewed as a fan-beam counterpart of the Katsevich algorithm for the conventional helical CT (Katsevich 2002a, 2002b, 2004). The second class (Pan et al 2004, Zhuang et al 2004, Noo et al 2004) may be viewed as a fan-beam counterpart of the Zou–Pan algorithms for the conventional helical CT (Zou and Pan 2004a, 2004b). The second aspect is to develop new methods of image reconstruction from partially truncated projection data (Clackdoyle et al 2004, Noo et al 2004). When the novel concept of a virtual scanning trajectory is introduced (Clackdoyle et al 2004), it has been demonstrated that the first class of super-short scan reconstruction algorithms (Noo et al 2002, Kudo et al 2002, Chen 2003) may be used to reconstruct images from a partially truncated projection data set. In this method, it seems that the data rebinning process is indispensable. This rebinning algorithm addresses the problem in which the projection data have been truncated only for some view angles. However, it has also been demonstrated that the partial data truncation problem may be solved by using the second class of the reconstruction algorithms (Noo et al 2004). In this case, it is not necessary to explicitly introduce the data rebinning process. A feature of this truncated scanning configuration is that projection data may or may not be truncated for each view angle. The necessary condition for exact reconstruction of an ROI is that the projection data corresponding to the rays passing through the ROI are untruncated in all views. In addition, a large class of new image reconstruction formulae for fan-beam projection data have also been discovered, and have been proven to be equivalent to Radon’s original formula (Clackdoyle and Noo 2004). In this paper, we consider a specific fan-beam scanning configuration with data truncation at every view angle: fan-beam projection data are acquired using an asymmetric detector that covers only one half of the image object (figure 1). We refer to this scenario as a constant view truncation (CVT) problem. This is a special case of the more general partial view data truncation problems that frequently occur in nuclear medicine (Chang et al 1995, Loncaric et al 1995, La Riviere et al 2001), radiation therapy (Hooper and Fallone 2002), and diagnostic x-ray CT (Hsieh et al 2004), where the detector is large enough to cover more than half of the image object but not large enough to cover the whole image object at some view angles. If the data are acquired from a complete circle, it is known that the data set is sufficient to reconstruct the whole image object (Chang et al 1995). However, the standard fan-beam filtered backprojection (FBP) image reconstruction algorithm does not provide a mathematically exact image reconstruction without an explicit data rebinning process from truncated fan-beam projections to non-truncated parallel-beam projections or from truncated fan-beam projections to non-truncated fan-beam projections. This new algorithm provides a means for exact reconstruction of the entire image object, without invoking a rebinning operation, from truncated data. In this paper, we will demonstrate that the newly developed fan-beam image reconstruction algorithm (Zhuang et al 2004) via filtering the backprojection image of differentiated projection data (FBPD) provides a mathematically exact reconstruction algorithm for the aforementioned CVT problem. Using the FBPD algorithm, data rebinning is no longer
Exact fan-beam image reconstruction algorithm
Focal Spot Trajectory
1807
right side
left side
Region of Image Support
(Co
nst Half ant Det Vie etec w T tor run ca
tion
) Full Detector
(No Truncation)
Figure 1. The specific data truncation problem which will be addressed in this paper is that in which projection data are collected from a detector of only half the size required for a non-truncated data. The detector is positioned asymmetrically for every given view angle. Note that the definitions of the left-hand side of the detector and the right-hand side of the detector are taken from the source position looking towards the detector.
needed. On the basis of an appropriate choice of filtering lines, a local ROI may be exactly reconstructed using the truncated projection data acquired from a scanning path which is less than 2π in angular coverage. 2. Brief review of the fan-beam FBPD reconstruction method: without data truncation For simplicity, we consider image reconstruction from fan-beam projections on an arc source trajectory with radius r (figure 2). The source trajectory is parametrized as y(t) = r(cos t, sin t),
t ∈ [ti , tf ],
(2.1)
where t parametrizes the view angle and is measured counter-clockwise with respect to the x-axis, ti and tf are the view angles corresponding to the starting and ending points of the scanning path. The newly developed FBPD image reconstruction scheme (Zhuang et al 2004) can be briefly summarized as below. For a given image point x, suppose a straight line, l, x )) that passes through the image point x intersects the source trajectory at two points y(t a ( x )). According to the Tuy data sufficiency condition (Chen 2003, Tuy 1983), the and y(t b ( x ), tb ( x )] are sufficient projection data from the scanning path within the angular range [ta ( for the reconstruction of the image point x. Throughout this section, we assume that only x ), tb ( x )] be used. This amounts to a special the projection data within the angular range [ta (
1808
S Leng et al
Figure 2. Definition of each of the vectors and angles used in deriving the FBPD reconstruction formula.
‘one or zero’ weighting scheme among a large class of other potential weighting schemes (Zhuang et al 2004). On the basis of this specific choice of the weighting function, the filtering must be performed along the straight line l in the space of the backprojection image x ), tb ( x )], one may construct a of the differentiated data. Using the projection data from [ta ( backprojection image Q( x , x ), tb (x ) ∂ ∂ 1 − [gm (γ = γx , t) − gm (γ = γx + π, t)], (2.2) Q( x , x ) = dt | x − y (t)| ∂t ∂γ ta ( x) where gm (γ , t) is the measured fan-beam projection data at view angle t and fan angle γ , where γ is measured from the iso-ray using the counter-clockwise convention. The values of γx and βˆ x are determined by (figure 2) γx = φx − t − π,
βˆ x = (cos φx , sin φx ),
(2.3)
where x is an arbitrary point in the backprojection image space. After the formation of a backprojection image associated with the image point x, the final image reconstruction is obtained through a Hilbert transform along the filtering line l. Namely, Q( x , x )|l 1 f ( x ) = − 2 dxl , (2.4) 2π l xl − xl where xl and xl are the projection of vectors x and x onto the filtering line l. For clarity, one may choose a coordinate system such that the x-axis is parallel to the filtering line l. In this case the above filtering process may be explicitly written as +∞ Q(x, y; x , y) 1 dx . (2.5) f (x, y) = − 2 2π −∞ x − x Note that the range of integration extends to infinity. Thus, we need to calculate the backprojection image values over an infinite range. This is mathematically possible as demonstrated in the above formulae (2.2) and (2.3). However, it is not computationally tolerable. Fortunately, for a compactly supported image function f ( x ) the image values of f ( x ) on the nonzero portion of the filtering line may be exactly calculated using a finite range Hilbert transform (Tricomi 1957, Milkhlin 1957). If we assume that the image to be reconstructed is nonzero only within a circle of radius R0 , and the filtering operation is
Exact fan-beam image reconstruction algorithm
1809
performed along a horizontal line (l) as indicated by equation (2.5), then the following formula is suggested to reconstruct image points ( x ): √ +R 1 R 2 − x 2 1 Q( x , x ) + C , dx (2.6) f (x, y0 ) = − 2 √ 2π x − x R 2 − x 2 −R where x = (x, y0 ), x = (x , y0 ), R R0 and C is a constant which is filtering line dependent. x , x ) are Using this formula, only the data in the range x ∈ [−R, +R] of the function Q( needed in the image reconstruction. This finite Hilbert transform reduces the computational load significantly since we now only need to calculate the values of the backprojection image within a finite range. To determine the constant C in equation (2.6), one can use the fact that the image values should be zero when the image point is outside the region of support (Zhuang et al 2004, ¯ y0 ) is outside the image support, then the Noo et al 2004). Therefore, if a point xc = (x, constant C can be determined as √ +R R 2 − x 2 Q( x c , x ). dx (2.7) C=− −x ¯ x −R In practice, several points along the filtering lines, which were also outside the region of compact support, were chosen to calculate the constant C. The mean value of several calculated C values was used in (2.6) to determine the image values for points along the same filtering line. An alternative means to determine the constant C is to integrate the both sides of equation (2.6) along the filtering line l. Namely, integration over the variable x is calculated. The integral on the left-hand side gives nothing but the projection value along the filtering line and it is denoted as gm (γ0 , t0 ). The integral over the variable x in the first term of the right-hand side of equation (2.6) vanishes. The integral of the second term is given by −C/2π . Thus, the constant C is determined by the following equation (Tricomi 1957): C = −2πgm (γ0 , t0 ).
(2.8)
The advantage of this method to determine the constant C is that there is no need to calculate the image values for the points outside the image support. Therefore, we only need the backprojection image values in the range [−R0 , +R0 ]. In summary, we emphasize the following important feature of the FBPD method given by equations (2.2) and (2.6): if the value of the backprojection image Q( x , x ) is correct for all points along a filtering line then any point on that filtering line may be accurately reconstructed. Therefore, it is crucial to obtain correct values of the backprojection image Q( x , x ) for all the points along the filtering line l. The correct values of the backprojection image Q( x , x ) x ), tb ( x )] provides a projection value for any are achieved when: every view angle t ∈ [ta ( given point x on a chosen filtering segment. In the case of non-truncated data as discussed in Clackdoyle and Noo (2004), in order to obtain the correct value of the backprojection image Q( x , x ), there is a degree of freedom in the selection of the scanning path of the backprojection region for a given filtering line. Namely, one may backproject the differentiated data from either of two given arcs (figure 3). One may also selectively combine the projection values of the conjugate rays from two arcs in figure 3 to form the correct backprojection image. The use of compensation from conjugate rays will be elaborated upon below. 3. Fan-beam FBPD image reconstruction: constant view truncation (CVT) In the standard fan-beam filtered backprojection (FBP) reconstruction (Kak and Slaney 1987) algorithm and the super-short scan image reconstruction algorithms developed in Noo et al
1810
S Leng et al
ROI
ROI
Filtering Line (l)
(a)
Filtering Line (l)
(b)
Figure 3. Options of angular range for a backprojection image from non-truncated projection data. In order to achieve the correct backprojection image for all points on the filtering line given above, one may select the data from either the source trajectory highlighted in (a) or the trajectory highlighted in (b).
(2002), Kudo et al (2002) and Chen (2003), the filtering step is conducted over all of the projection data for each given view angle. The filtered data are then backprojected in order to form the reconstructed image. Note that the filtering kernels, either a ramp filter or a Hilbert filter, in the aforementioned FBP-type image reconstruction algorithms are non-local. Therefore, for the truncated projection data, an explicit rebinning process is necessary for the use of the previously mentioned FBP image reconstruction algorithms. Originally, the FBPD image reconstruction algorithm (Zhuang et al 2004) was derived under the assumption that the projection data are not truncated for any view angle. In this section, we will elaborate on this point carefully to demonstrate that the final FBPD image reconstruction algorithm given by equations (2.2) and (2.6) survives the CVT problem. In order to address the CVT problem using the FBPD algorithm, the following observations are crucial: (1) for a given point x in the backprojection image space, only one projection for each view angle is required to construct the backprojection image. Therefore, for a given x ), tb ( x )] and a given point x in the backprojection image space, one view angle t ∈ [ta ( does not have to know all of the projection values within the fan-shape region that covers the whole image object in order to properly reconstruct the image. (2) For a given point x in the backprojection image space, only one ray is required in (1). This ray may be from a view x ), tb ( x )] or it may be a conjugate ray from a view angle out of the angular range angle t ∈ [ta ( x ), tb ( x )]. The second observation implies that one may combine the projection data from [ta ( different view angles in order to form a correct backprojection image. We analyse below how to apply the FBPD reconstruction algorithm to the CVT problem (figure 1), where an asymmetrically placed half-size detector is utilized to acquire fan-beam projection data. In order to facilitate the analysis, we assume that the image function has a definite boundary (figure 4). Outside the boundary , the image function f ( x ) ≡ 0. For a given image point x, the filtering line l intersects the boundary curve at a minimum of two points A and B (figure 4). The detector is large enough to cover one half of the image object. Namely, if the iso-ray divides the image object into two parts at any view angle, the detector is always sufficiently large to guarantee that there are no data truncation for one of the two parts. Each of the schematic diagrams in the present paper assumes the convention that the data provided are on the right-hand side of the iso-ray given that one is oriented from the source towards the detector (figure 1). A third-generation scanning configuration and an equal-angular detector
Exact fan-beam image reconstruction algorithm
l A
x
B
1811
A
B
l
x
Γ
Γ
(a)
(b)
Figure 4. Examples of possible filtering lines l, and the definition of points A and B which are the outer most points intersection points between the filtering line and the curve . Source Trajectory Illuminating Sample Image Point
F2
F1
F3
x x
(a)
x
x
(b)
(c)
(d)
Figure 5. Demonstration that in the CVT case any one given image point is only illuminated by rays originating from an arc of arc-length π . Three sample focal spots F1 , F2 and F3 are given to illustrate this point in (a), (b) and (c); the total angular range is shown in (d).
are also assumed. The analysis may easily be adapted to the case of a collinear detector or the case of a left-hand side convention of object illumination. Using the above scanning configuration, for a given image point x, it is easy to see that only one half of the circular scanning path will illuminate the image point x (figure 5). Geometrically, if one draws a straight line that passes through the points x and the iso-centre, the straight line bisects the circular scanning path. According to the convention of the scanning configuration, one may easily determine which half of the scanning path will illuminate the point x. We are now in the position to discuss the impact of the selection of the filtering line on the reconstruction algorithm. The choice of filtering lines in the reconstruction algorithm will determine the angular range of data required to achieve an exact reconstruction for all of the points on that given line. In order to provide a general foundation for the selection of filtering lines two distinct classes of filtering lines will be discussed below. The first class is all the filtering lines which pass through the origin, and the second class is the filtering lines which do not pass through the origin. The angular range of data required for each of these scanning paths will also be addressed. Following this analysis will be a discussion of how one may reconstruct an ROI from a group of filtering lines. In the first class of filtering lines, one selects a filtering line as the line that connects an image point x and the iso-centre (figure 4(a)). Note that one half of the scanning path will illuminate boundary point A. In order to illuminate the other boundary point B, data are required from the other half of the scanning path. Therefore, only a full scan can generate
1812
S Leng et al Required Source Trajectory
ROI
Required Source Trajectory ROI
Correct Backprojection Data
ROI
Filtering Line (l)
Correct Backprojection Data
Correct Backprojection Data
Required Source Trajectory
(a)
(b)
(c)
Figure 6. In comparison with the case of untruncated data in figure 3 in order to provide the correct backprojection data for a central filtering line a full scan is needed for the case of constant view truncation. In this case it is evident that the union of the backprojection data required in (a) and (b) form the complete 2π trajectory (c).
sufficient projection data for all the points along the radial line (figure 6). This point may be further clarified by analysing the illumination pattern given in figures 5(a)–(c). When the focal spot moves from figure 5(a) to figure 5(c), only the left-hand side of the horizontal line is illuminated. In other words, one may obtain the correct values of the backprojection image for all the points on the left portion of the horizontal filtering line (figures 5(d) and 6(a)). In order to obtain the correct backprojection image values required for the right portion of the horizontal filtering line, one must scan through the other half of the circular scanning path (figure 6(b)). Therefore, a full 2π angular range is required to provide sufficient data to generate the correct backprojection images for all of the points lying along this filtering line (figure 6(c)). This situation is significantly different from the non-truncated case illustrated in figure 3. Due to the symmetry breaking when moving from the non-truncated data reconstruction problem to the half-size detector CVT reconstruction problem, the degree of freedom in the backprojection step to choose either arc (figure 3(a) or figure 3(b)) has been lost. In the CVT reconstruction case in order to generate correct values for the backprojection image at each point along this radial filtering line data are required from the full angular range (figure 6(c)). The second class of filtering lines is composed of those filtering lines that do not pass through the iso-centre. One such example is illustrated in figure 7. In this case, two half circle illumination regions corresponding to the boundary points A and B overlap with one another. The union of these illumination regions is less than a full circle. In the radial filtering case discussed above the union of these two illumination regions consisted of the full 2π trajectory (figure 6). The starting and ending view angles will be denoted as ti and tf . It is important to observe that the view angles in the angular range [ta , ti ], which are formally required by the FBPD reconstruction algorithm, are missing. Thus, there are no rays emanating from this angular range. The data in the angular range [ta , ti ], which are required in order to generate the correct values of the backprojection image associated with point A are not present in the region defined by the union of the illumination regions for points A and B. However, one can use the conjugate rays emanating from the view angles in the angular range [ti , tf ] to compensate for the missing projections in the angular range [ta , ti ] (figure 8). Therefore, the angular range [ti , tf ] is sufficient for the reconstruction of the filtering segment AB . To briefly summarize the analysis given above, we conclude that the selection of the filtering lines is an important issue and has implications on the required angular range of projection data. Now that the effect of the orientation of the filtering lines has been probed we may address the issue of reconstructing an ROI from a combination of filtering lines. In the case of radial
Exact fan-beam image reconstruction algorithm
1813
l
A
B
x
Γ
Figure 7. In order to reconstruct x and all other points along the filtering line l the required source trajectory is given by the union of the all source positions that illuminate image point A and all the source positions that illuminate image point B. Conjugate Ray Conjugate Source Position
l
Missing Source Position
ti ta
Conjugate Rays
B
A
Γ
l
ti ta
B
A
Γ
tf
tf (a)
(b)
Figure 8. The data in the angular range [ta , ti ] are missing. However, since for each view angle the backprojection image only has a contribution from a single projection, the missing projections may be compensated for by substituting their equivalent conjugate projections. The conjugate rays are emitted from different view angles (a). In each case the conjugate source point is within the union of the the illumination regions of points A and B (b). Remember that the right-hand convention is still used here and thus it is obvious that the conjugate ray will not be truncated.
filtering (figure 4(a)), there is no small ROI reconstruction. Since the union of the regions of illumination for any given radial filtering lines composes the entire 2π circle (figure 6), there is no possibility of a ‘short scan’. However, in the case of non-radial filtering lines (figure 4(b)) the possibility of a ‘short scan’ exits. As shown in figure 9 the angular range required for the outermost filtering line of the scanning path is sufficient to provide a complete backprojection data set for an ROI bounded by the filtering line and the boundary curve . Thus, for non-radial filtering lines one may accurately reconstruct a local ROI using a ‘short scan’ mode. Namely, the angular range of
1814
S Leng et al Required Source Trajectory
Reconstructed ROI B
l
B
B
l
A
B
l l
A
A
Γ
(a)
Γ
(b)
Γ
(c)
A
Γ
(d)
Figure 9. Demonstration that the angular range required to generate the backprojection data for the outermost filtering line (a) is sufficient to generate the backprojection data for an ROI. This illustrative case shows that an ROI may be reconstructed using filtering lines parallel to the outermost filtering line (b) and (c). Thus, the FBPD provides the opportunity for a ‘short scan’ from CVT data (d).
the scanning path is less than 2π . Note that the meaning of our short scan is different from the conventional case with scanning angular range π + γm (γm is the fan-angle of the whole image object) (Kak and Slaney 1987). The selection of the optimal filtering lines is dependent upon the given scanning configuration and the desired ROI, and thus the topic of the selection of the optimal filtering lines will not be addressed here. It is important to emphasize that the short-scan mode is a significant result of the FBPD based image reconstruction algorithm. In the conventional strategy of rebinning the truncated fan-beam projections to non-truncated parallel beam projections, no short-scan mode is possible. Namely, even for the case of small ROI reconstruction, data from the entire scanning path are required. 4. Numerical simulations In this section, we present numerical simulations in order to illustrate the validity of the proposed fan-beam FBPD image reconstruction from the constantly truncated projection data. A standard Shepp–Logan phantom (Kak and Slaney 1987) was used to generate the projection data from an arc source trajectory. The Shepp–Logan dimensions were scaled by a constant factor of 5.8. The total fan angle was chosen to be γm = π4 . The sampling rate of the equiγm . The sampling rate of the view angles was chosen angular detector was chosen as γ = 723 π to be t = 1031 . Each of the images displayed in this section has a compressed window of [0.98 1.06] in order to highlight the low contrast objects in the phantom. In computing the backprojection image of differentiated projection data, the derivatives were calculated by a three-point formula. A Hann window was used to suppress the highfrequency fluctuations in the reconstructed images. The constant C in equation (2.6) was determined by using equation (2.8). A bilinear interpolation scheme was used to find the projection values for a given filtering line. Figure 10 shows a simulation for a complete scan, where horizontal parallel filtering lines were chosen. Figure 10(a) demonstrates the scanning path and orientation of the filtering lines. Figure 10(b) shows the reconstructed image. Figures 10(c) and 10(d) present the central horizontal line and the central vertical line of the image in figure 10(b). Theoretical values are presented for comparison in each case. The simulation results show that the whole object can be reconstructed correctly when a full scan is used with the asymmetric half detector.
Exact fan-beam image reconstruction algorithm
A
1815
B
Full FOV Supported (a) 1.08
(b)
Theoretical Values Reconstructed Values
Theoretical Values Reconstructed Values
1.06
1.06 1.05 1.04 1.04 1.02 1.03 1
1.02
0.98
1.01
0.96
1 –4
–3
–2
–1
0
(c)
1
2
3
4
–6
–4
–2
0
2
4
6
(d)
Figure 10. Simulation results of a full scan (2π ) source trajectory. A schematic figure illustrating the scanning path used and the orientation of filtering lines (a). Parallel horizontal filtering lines were used here. The reconstructed image (b) along with density plots of the central horizontal and vertical lines (c) and (d) demonstrate good agreement with theoretically expected values.
Figures 11 and 12 are simulation results of reconstructing a specific ROI using a short scan. The scanning path was parametrized as (2π − α), where α = π/4 for each of the ROI simulations given below. Figures 11(a) and 12(a) demonstrate the scanning path and direction of the filtering lines for ROIs in the upper and lower portion of the image respectively. Figures 11(b) and 12(b) demonstrate the reconstructed images. Figures 11(c) and 12(c) present density plots from horizontal lines inside of ROI, and 11(d) and 12(d) present horizontal lines outside of ROI. The density values within the ROI have been accurately reconstructed. Quantitatively exact CT number cannot be obtained outside the ROI, but some of the image structure remains visible in regions near the ROI. Figures 11(e) and 12(e) plot the density values along central vertical lines of the reconstructed images. The region of support for exact reconstruction is defined by (y > 1.6) and (y < −1.6) for figures 11 and 12 respectively. In each case the theoretical phantom values are presented in each plot for comparison. Figure 13 is the result of a short scan where the ROI is the left most portion of the phantom. The main difference between this reconstruction and that presented in
1816
S Leng et al
Supported ROI B
A
γm
(a)
(b) 1.08
Theoretical Values Reconstructed Values
1.08
Theoretical Values Reconstructed Values
1.06 1.06 1.04 1.04 1.02
1.02
1
1
0.98
0.98
0.96
0.96 –4
–3
–2
–1
0
1
2
3
4
–4
(c)
–3
–2
–1
0
1
2
3
4
(d)
Theoretical Values Reconstructed Values
1.06
1.05
1.04
1.03
1.02
1.01
1
–6
–4
–2
0
2
4
6
(e)
Figure 11. Reconstruction of an ROI from the top part of the object using a short-scan source trajectory. A schematic figure of this short-scan trajectory and orientation of the outermost filtering line is given in (a). Parallel horizontal filtering lines were used to reconstruct the ROI. The dashed arc presents the angular range of source trajectory. Reconstructed image where the ROI is given by the region above the solid line (b). Density plot of two horizontal lines (c) and (d) at the positions indicated by the dashed lines in (b), where (c) is within the ROI and (d) is outside the ROI. The density plot of the central vertical line is given in (e), where the image values for the exact ROI correspond to y > 1.6.
Exact fan-beam image reconstruction algorithm
1817
γm A
B
Supported ROI (a)
(b)
Theoretical Values Reconstructed Values
1.08
1.08
1.06
1.06
1.04
1.04
1.02
1.02
1
1
0.98
0.98
Theoretical Values Reconstructed Values
0.96
0.96 –4
–3
–2
–1
0
1
2
3
4
(c) 1.06
–4
–3
–2
–1
0
1
2
3
4
(d)
Theoretical Values Reconstructed Values
1.05
1.04
1.03
1.02
1.01
1
–6
–4
–2
0
2
4
6
(e)
Figure 12. Reconstruction of an ROI from the bottom portion of the phantom using a short-scan source trajectory. A schematic figure of this short-scan trajectory and orientation of the outermost filtering line is given in (a). Parallel horizontal filtering lines were used to reconstruct the ROI. The dashed arc presents the angular range of source trajectory. Reconstructed image where the ROI is given by the region below the solid line (b). Density plot of two horizontal lines (c) and (d) at the positions indicated by the dashed lines in (b), where (c) is within the ROI and (d) is outside the ROI. The density plot of the central vertical line is given in (e), where the image values for the exact ROI correspond to y < −1.6.
1818
S Leng et al
A
Supported ROI
B
γm
(a) 1.08
(b)
Theoretical Values Reconstructed Values
Theoretical Values Reconstructed Values
1.08
1.06 1.06 1.04 1.04 1.02 1.02 1 1 0.98 0.98 0.96
0.96 –4
–3
–2
–1
0
1
2
3
4
(c)
1.08
–5
–4
–3
–2
–1
0
1
2
3
4
5
(d)
Theoretical Values Reconstructed Values
1.06
1.04
1.02
1
0.98
0.96 –4
–3
–2
–1
0
1
2
3
4
(e)
Figure 13. Reconstruction of an ROI from the left portion of the object using a short-scan source trajectory. A schematic figure of this short-scan trajectory and orientation of the outermost filtering line is given in (a). Parallel vertical filtering lines were used to reconstruct the ROI. The dashed arc presents the angular range of source trajectory. Reconstructed image where the ROI is given by the region to the left of the solid line (b). Density plot of two vertical lines (c) and (d) at the positions indicated by the dashed lines in (b), where (c) is within the ROI and (d) is outside the ROI. The density plot of the central horizontal line is given in (e), where the image values for the exact ROI correspond to x < −1.6.
Exact fan-beam image reconstruction algorithm
1819
figures 11 and 12 is that in this case vertical filtering lines were chosen. Figure 13(a) demonstrates the scanning path and orientation of the filtering lines. Figure 13(b) demonstrates the reconstructed image. Figures 13(c) and (d) present vertical lines which lie both inside and outside of the ROI. Figure 13(e) presents the density values of the central horizontal line in the reconstructed image. Note that in figures 11(e), 12(e) and 13(e) as expected the reconstructed values match the theoretical values within each ROI, and as the distance increases from the ROI the reconstructed values become increasingly less accurate. In summary fan-beam image reconstruction simulations were performed using the FBPD reconstruction algorithm on CVT data for scans of both the entire object as well as specific ROI reconstructions. In the case of ROI reconstructions a ‘short scan’ was used in which the angular range was less than 2π . The simulations support the theoretical assertions made above, and demonstrate that the FBPD algorithm provides accurate reconstruction of the entire object (full-scan case) or of a given region of interest (short-scan case). 5. Conclusions and discussion In this work, we present a new image reconstruction algorithm for a specific data truncation problem. The scanning configuration is that of a third generation geometry CT scanner in which only a detector of half of the size required for non-truncated projections is placed asymmetrically about the iso-ray. In this geometry, the data are truncated at every view angle, and thus this problem has been referred to as constant view truncation. If the scanning path is a complete circle, it was previously demonstrated that the data are sufficient to reconstruct the whole image object (Chang et al 1995) after the fan-beam projection data are rebinned into parallel-beam projection data. In this rebinning scheme, even in the case of ROI reconstruction a scanning path shorter than 2π is not sufficient. In this work, we demonstrate that an ROI image may be accurately reconstructed using constantly truncated projection data collected from an angular range of less than the full 2π scan. The connection between the acquired data and the generated backprojection image, in the recently developed FBPD reconstruction (Zhuang et al 2004), has been exploited to yield the feasibility of a ‘short scan’ in the constant view truncation (CVT) reconstruction problem. Numerical simulations were conducted to validate the proposed image reconstruction scheme for both an ROI reconstruction and the reconstruction of a whole object. Acknowledgments This work was partially supported by National Institute of Health grants 1R21 EB001683-01, 1R21 CA109992-01, 5T32CA009206-24, a start-up grant from the University of Wisconsin– Madison, and a grant from GE Healthcare Technologies. References Chang W, Loncaric S, Huang G and Sanpitak P 1995 Asymmetric fan transmission CT on SPECT system Phys. Med. Biol. 40 913–28 Chen G-H 2003 A new framework of image reconstruction from fan beam projections Med. Phys. 30 1151–61 Clackdoyle R and Noo F 2004 A large class of inversion formulae for the 2D Radon transform of functions of compact support Inverse Problems 20 1281–91 Clackdoyle R, Noo F, Guo J and Roberts J 2004 A quantitative reconstruction from truncated projections in classical tomography IEEE Trans. Nucl. Sci. 51 2570–8 Hooper H R and Fallone B G 2002 Sinogram merging to compensate for truncation of projection data in tomotherapy imaing Med. Phys. 29 2548–51
1820
S Leng et al
Hsieh J, Chao E, Thibault J, Grekowicz B, Horst A, McOlash S and Myers T J 2004 A novel reconstruction algorithm to extend the CT scan field-of-view Med. Phys. 31 2385–91 Kak A C and Slaney M 1987 Principles of Computerized Tomographic Imaging (New York: IEEE Press) Katsevich A 2002a Theoretically exact filtered backprojection-type inversion algorithms for spiral CT SIAM J. Appl. Math. 62 2012–26 Katsevich A 2002b Analysis of an exact inversion algorithm for spiral cone-beam CT Phys. Med. Biol. 47 2583–97 Katsevich A 2004 An improved exact filtered backprojection algorithm for spiral computed tomography Adv. Appl. Math. J. Math. Math. Sci. 32 681–97 Kudo H, Noo F, Defrise M and Clackdole R 2002 New super-short scan algorithms for fan-beam and cone-beam reconstruction Conf. record of the 2002 IEEE Nuclear Science Symposium and Medical Imaging Conference (Norfolk, VA) (IEEE Service Centre vol 2) pp 902–6 La Riviere P J, Pan X, Gilland D, Kao C-M, Chang W and Jaszczak R 2001 Transmission image reconstruction and redundant information in SPECT with asymmetric fanbeam collimation IEEE Trans. Nucl. Sci. 48 1357–63 Loncaric S, Chang W and Huang G 1995 A processing technique for the truncated projections of asymmetric-fan-beam transmission imaging IEEE Trans. Nuc. Sci. 42 2292–7 Milkhlin S G 1957 Integral Equations (New York: Pergamon) pp 126–31 Noo F, Defrise M, Clackdoyle R and Kudo H 2002 Image reconstruction from fan-beam projections on less than a short-scan Phys. Med. Biol. 47 2525–46 Noo F, Clackdoyle R and Pack J D 2004 A two-step Hilbert transform method for 2D image recontruction Phys. Med. Biol. 49 3903–23 Pan X, Xia D, Zou Y and Yu L 2004 A unified analysis of FBP-based algorithms in helical cone-beam and circular cone- and fan-beam scans Phys. Med. Biol. 49 4349–69 Tricomi F G 1957 Integral Equations (New York: Interscience) pp 173–85 Tuy H K 1983 An inversion formula for cone-beam reconstruction SIAM J. Appl. Math. 43 546–52 Zhuang T, Leng S, Nett B E and Chen G-H 2004 Fan-beam and cone-beam image reconstruction via filtering the backprojection image of differentiated projection data Phys. Med. Biol. 49 5489–503 Zou Y and Pan X 2004a Exact image reconstruction on PI-lines from minimum data in helical cone-beam CT Phys. Med. Biol. 49 941–59 Zou Y and Pan X 2004b Image reconstruction on PI-lines by use of filtered backprojection in helical cone-beam CT Phys. Med. Biol. 49 2717–31
