Use of Hilbert transforms to simplify image reconstruction in a spiral scan cone beam CT imaging system

[0002] 1. Field of the Invention

[0003] The invention relates generally to a cone beam computed tomography (CT) imaging system, and more specifically to a simplified method and apparatus for image reconstruction in such a system.

[0004] 2. Description of the Prior Art

[0005] Recently a system employing cone beam geometry has been developed for three-dimensional (3D) computed tomography (CT) imaging that includes a cone beam x-ray source and a 2D area detector. An object to be imaged is scanned, preferably over a 360° angular range and along its entire length, by any one of various methods wherein the position of the area detector is fixed relative to the source, and relative rotational and translational movement between the source and object provides the scanning (irradiation of the object by radiation energy). The cone beam approach for 3D CT has the potential to achieve 3D imaging in both medical and industrial applications with improved speed, as well as improved dose utilization when compared with conventional 3D CT apparatus (i.e., a stack of slices approach obtained using parallel or fan beam x-rays).

[0006] As a result of the relative movement of the cone beam source to a plurality of source positions (i.e., “views”) along the scan path, the detector acquires a corresponding plurality of sequential sets of cone beam projection data (also referred to herein as cone beam data or projection data), each set of cone beam data being representative of x-ray attenuation caused by the object at a respective one of the source positions.

[0007] As well known, and fully described for example in the present inventor's U.S. Pat. No. 5,257,183 entitled METHOD AND APPARATUS FOR CONVERTING CONE BEAM X-RAY PROJECTION DATA TO PLANAR INTEGRAL AND RECONSTRUCTING A THREE-DIMENSIONAL COMPUTERIZED TOMOGRAPHY (CT) IMAGE OF AN OBJECT issued Oct. 26, 1993, incorporated herein by reference, image reconstruction processing generally begins by calculating Radon derivative data from the acquired cone beam data. The Radon derivative data is typically determined by calculating line integrals for a plurality of line segments L drawn in the acquired cone beam data. In the embodiment described in detail in the 5,257,183 patent, Radon space driven conversion of the derivative data is used to develop an exact image reconstruction of a region-of-interest (ROI) in the object.

[0008] A cone beam data masking technique which improves the efficiency of the calculation of the Radon derivative data in such a Radon space driven technique is described in the present inventor's U.S. Pat. No. 5,504,792 entitled METHOD AND SYSTEM FOR MASKING CONE BEAM PROJECTION DATA GENERATED FROM EITHER A REGION OF INTEREST HELICAL SCAN OR A HELICAL SCAN, issued Apr. 2, 1996, also incorporated herein by reference. The masking technique facilitates efficient 3D CT imaging when only the ROI in the object is to be imaged, as is normally the case. In the preferred embodiment described therein, a scanning trajectory is provided about the object, the trajectory including first and second scanning circles positioned proximate the top and bottom edges, respectively, of the ROI, and a spiral scanning path is connected therebetween. The scanning trajectory is then sampled at a plurality of source positions where cone beam energy is emitted toward the ROI. After passing through the ROI the residual energy at each of the source positions is acquired on an area detector as a given one of a plurality of sets of cone beam data. Each set of the cone beam data is then masked so as to remove a portion of the cone beam data that is outside a given sub-section of a projection of the ROI in the object and to retain cone beam projection data that is within the given sub-section. The shape of each mask for a given set of cone beam data is determined by a projection onto the detector of the scan path which is above and below the source position which acquired the given set of cone beam data. The masked (i.e., retained) cone beam data is then processed so as to develop reconstruction data An exact image of the ROI is developed by combining the reconstruction data from the various source positions which intersect a common integration plane. Hence, the masks are commonly referred to as “data-combination” masks.

[0009] Data-combination masks can also be used to improve the efficiency of the calculation of the derivative data in a detector data driven technique, such as the simplified ramp filter technique described in the present inventor's U.S. Pat. No. 5,881,123 entitled SIMPLIFIED CONE BEAM IMAGE RECONSTRUCTION USING 3D BACKPROJECTION, issued Mar. 9, 1999, also incorporated herein by reference. This simplified technique reconstructs the image using 2D approximation data sets formed by ramp filtering the masked cone beam data in the direction of the projection of a line drawn tangent to the scan path at the source position that acquired that set of cone beam data. Although this technique is less complex than the prior techniques, the reconstructed image is not exact.

[0010] Accordingly, U.S. patent application Ser. No. 09/162,303 entitled ADAPTIVE MASK BOUNDARY CORRECTION IN A CONE BEAM

[0011] IMAGING SYSTEM, filed Sep. 28, 1998, now U.S. Pat. No. 6,084,937 issued Jul. 4, 2000, and also incorporated herein by reference, describes a technique for computing 2D correction data which, when combined with the ramp filtered 2D approximation data sets, yields an exact image reconstruction. The 2D correction data basically comprises a point spread function representative of image reconstruction processing for each point on the detector which interesects the boundary of the data-combination mask. However, the techniques disclosed therein require circle scan paths near the top and bottom edges of the ROI when only an ROI is to be imaged. It would be desirable to provide such a simplified image reconstruction processing without the requirement of circle scan paths near the top and bottom edges of the ROI.

[0012] The present inventor's allowed U.S. patent application Ser. No. 09/343,770 entitled EXACT REGION OF INTEREST CONE BEAM IMAGING WITHOUT CIRCLE SCANS, filed Jun. 30, 1999, incorporated by reference herein, improved upon the invention described in the forenoted U.S. Pat. No. 5,504,792, by providing an exact image reconstruction of an ROI in an object without the requirement that the source scan path have top and bottom circle scan paths near the top and bottom edges of the ROI. Furthermore, the improvement is applicable to both of the Radon space and detector driven types of image reconstruction processing. As described in this U.S. patent application Ser. No. 09/343,770, and consistent with the techniques described in the above noted U.S. Pat. Nos. 5,881,123 and 5,504,792, when calculating the derivative data, the length of the line segments L formed in the acquired cone beam data (along which integral data is to be developed) are determined by the boundaries of the data-combination mask. Exact image reconstruction processing for source positions which are “internal” to the top and bottom edges of the ROI can use the simplified ramp filtering with 2D correction data technique of the forenoted U.S. patent application Ser. No. 09/162,303. However, when processing cone beam data acquired at source positions near the top or bottom edges of the ROI, the line segments L are divided into two groups, where in one of the groups the line segments L have one of their end points determined by a horizontal line of the mask (in a preferred embodiment this horizontal line is the x-axis of the mask), instead of having both of the endpoints determined by the outer boundaries of the mask. When processing along the line segments L in this one group, acquired cone beam data which resides along lines L on the other side of the horizontal axis of the mask are not used. Thereafter, integral data is calculated for the line segments L having their endpoints limited by the horizontal axis in the mask, and the data is further processed so as to develop contributions to a 3D image reconstruction of the ROI in the object.

[0013] More specifically, when processing cone beam data acquired at source positions near the top or bottom edges of the ROI, the Filtered Backprojection (FBP) image reconstruction technique described in the forenoted U.S. patent application Ser. No. 09/343,770 can be used. FBP image reconstruction consists of two different kinds of processing: the first kind is 2-dimensional (2D) and the second kind is 3-dimensional (3D). In the 2D step each cone beam projection image is processed in a 2D space, and in the 3D step the processed 2D images are backprojected into a 3D object space. As shown in FIGS. 1A and 1B, the 2D step consists of the following 4 FBP image reconstruction sub-steps for processing the cone beam data acquired at each of a plurality of the source positions (S_i) along the scan path:

[0014] 1. Compute a 1D projection (i.e., a line integral) of the cone beam image acquired on a detector plane 100, at each of a plurality of angles θ. This step is illustrated in FIG. 1A for a given angle θ₁of a plurality of angles θ₁. A 1D projection 102 is shown at coordinates r, θ₁comprising the integrated values of the cone beam image 104 on detector plane 100 along a plurality of parallel lines L(r, θ₁) that are normal to angle θ₁, each line L being at an incremental distance r from an origin O. As shown and described below in conjunction with FIGS. 3A and 3B illustrating processing of cone beam data acquired near the top edge of the ROI, the lengths of the lines L will be limited using the masking techniques of the forenoted U.S. patent application Ser. No. 091343,770. Generally, if the detector plane 100 comprises an N by N array of pixels, then the number of angles θ is typically given by πN/2.

[0015] 2. Filter (differentiate) each 1D projection 102 in accordance with a d/dr filter, resulting in a new set of values at each of the r, θ coordinates, such as shown by the derivative projection 106 for the angle θ₁, shown in FIG. 1A. Note, the sum (integration) of the resulting values at these r, θ coordinates yield a quantity proportional to the Radon derivative for an integration plane Q(r, θ), as described above for Radon space driven image reconstruction processing.

[0016] 3. As illustrated by FIG. 1B, backproject the derivative projection 106 from each angle θ into a 2D object space 107 (which coincides with the detector plane 100). Lines 108 are representative of this backprojection, and spread the value from each r coordinate into 2D space 107 in a direction normal to each θ, thereby developing contributions to a backprojection image 109. Note, 2D object space 107 has a size corresponding to a virtual detector which is enlarged (compared with detector having a height corresponding to the pitch of the scan path), so as to cover the entire ROI in the object. This enlargement is required because the calculated Radon data affects the reconstruction of the entire Q plane, and not just the partial plane represented by the data combination mask.

[0017] 4. Perform a 1D d/dt filtering of the backprojection image formed in 2D space 107 by step 3. The 1D filtering is performed in the direction of the scan path, i.e., along lines 110, where t points in the direction of the projection of a line drawn tangent to the scan path.

[0018] As shown in FIG. 1B, the 3D step comprises performing a weighted 3D backprojection of the resulting data from 2D space 107 (i.e., from each pixel in the detector) onto a plurality of sample points P in a 3D object volume 112. The density assigned to each point P is weighted by the inverse of the square of the distance between the sample point and the location of the x-ray.

[0019] Although this technique allows exact image reconstruction without the requirement of a circle scan path near the top and bottom edges of the ROI in the object, the 4-substeps of 2D processing as described above, is relatively complex and time consuming, since it is required to be repeated for each angle θ.

[0020] It would be desirable to provide a simplified (i.e., less complex and time consuming) image reconstruction processing, such as that described in the forenoted U.S. Pat. No. 5,881,123, which does not require the use of circle scans near the top and bottom edges of the object being imaged.

[0021] A scanning and data acquisition method and apparatus for three dimensional (3D) computerized tomography (CT) imaging of an ROI in an object, wherein a reconstructed image is develop by calculating reconstruction data along a plurality of line segments L formed in cone beam data acquired on a detector at a plurality of source positions about the ROI. The endpoints of the line segments L formed in the cone beam data acquired at each of the source positions is determined by a data combination mask which is applied to the cone beam data. At source positions near the top and bottom of the ROI, the mask is divided into separate portions using at least a horizontal line therein, each portion having a corresponding set of line segments L formed therein. The reconstruction data is calculated by performing a plurality of one-dimensional (1D) convolving operations of the masked projection data along a respective set of the line segments L, the 1D convolving operation for each set of the line segments L developing an additive contribution to a new projection image. Finally, the new projection image is backprojected into a 3D space, thereby reconstructing a 3D image of the ROI in the object.

[0022] Thus, in accordance with the principles of the present invention, the most time-consuming tasks in backprojection ROI spiral scan cone beam CT, i.e., the 2D projection and 2D backprojection on the detector plane for the views near the two ends of the spiral, are combined into 1D Hibert Transforms, thereby tremendously simplifying the image reconstruction operation.

[0023] FIGS. 1A and 1B, previously described, illustrate a prior art 3D backprojection approach for reconstructing an image using the apparatus of FIG. 2;

[0024] FIG. 2 illustrates a 3D CT imaging apparatus useful for acquiring and processing cone beam data in accordance with the principles of the present invention;

[0025] FIGS. 3A and 3B illustrate masking of lines formed in the sets of cone beam data acquired near the top edge of an ROI in an object, for use reconstructing an image in accordance with a prior art technique;

[0026] FIG. 4 illustrates a top view of a portion of the image reconstruction geometry used to determine the point C₀in the masks shown in FIGS. 3, 7 and 9 to 11;

[0027] FIG. 5 illustrates geometry usefull for understanding how the combined operations of FBP image reconstruction sub-steps (1)-(3) of FIGS. 1A and 1B are equivalent to a 1D Hibert transform in the direction t;

[0028] FIGS. 6A, 6B and 6C illustrate a geometry useful for understanding how the combined operations of FBP image reconstruction sub-steps (1)-(3) within a given angular range are equivalent to the sum of two 1D Hilbert transforms.

[0029] FIG. 7 illustrates a portion of image reconstruction processing in accordance with the invention for a source point P which is below the x-axis in the data-combination mask;

[0030] FIGS. 8A and 8B illustrate geometry useful for understanding how the image reconstruction processing varies depending upon where the source point P resides on the top half of the data-combination mask;

[0031] FIGS. 9A, 9B and 9C illustrate portions of image reconstruction processing in accordance with the invention for a source point P which is positioned in the data-combination mask according to the criteria of FIG. 8A;

[0032] FIGS. 10A, and 10B illustrate portions of image reconstruction processing in accordance with the invention for a source point P which is positioned in the data-combination mask according to the criteria of FIG. 8B;

[0033] FIG. 2 illustrates a cone beam 3D CT imaging apparatus useful for acquiring and processing acquired cone beam data in accordance with the principles of the present invention. The illustrated imaging apparatus is constructed and operates substantially in accordance with the same principles described in the forenoted U.S. Pat. No. 5,257,183 and uses a masking image reconstruction technique such as described in the forenoted U.S. patent application Ser. No. 09/343,770, already incorporated herein by reference. Furthermore, image reconstruction processing in accordance with the present invention is accomplished in conjunction with a radiation blocking shutter for blocking radiation exposure to a portion of the object at source positions near the top and bottom edges of the ROI, as will be more specifically described below.

[0034] As shown in FIG. 2, a computer controlled manipulator 8, in response to control signals from an appropriately programmed computer 6, cause a source 10 of a cone or pyramid shaped beam of energy (such as x-rays) and a two-dimensional pixelated detector array 12 to cooperate (scan) at a plurality of discreet, sequentially occurring adjacent source positions, along a pre-defined source scanning path. In the illustrated embodiment the scanning path is shown as a spiral scan path 14 centered on a predetermined axis 15 of an object 16, and which proceeds in a counter-clockwise direction and spiraling upward. Other types of scan paths that encircle and traverse object 16 can also be used, however, a scan path 14 exhibiting a high degree of symmetry in its parallel projection is preferred.

[0035] The only height requirement on the detector is that it should be more than the distance between adjacent turns of a cone beam projection of the spiral scan path 14 onto the detector 12. As described in the forenoted U.S. patent application Ser. No. 09/343,770, when only a “region of interest” (ROI, shown as the shaded area of object 16) is to be imaged, providing top and bottom scan path circles at the top and bottom ends, respectively, of the ROI, are not necessary. Thus, a simple continuous spiral scan path 14 that traverses the length of the ROI is all that is needed.

[0036] As a result of the source/detector cooperation under control of computer 6 and manipulator 8, at each of the source positions, beams of x-ray energy pass through the field of view of the imaging apparatus, are selectively attenuated by varying energy absorbing densities inside object 16, and a set of cone beam data corresponding to the sensed x-ray energy falling on the elements (pixels) of detector 12 are developed. The sets of cone beam data are then supplied to a data acquisition system (DAS) 17 which, like the previously described portions of FIG. 1, may operate in a fashion well known to those of ordinary skill in this technology for digitizing and storing of the acquired cone beam data.

[0037] As described in the present inventor's forenoted U.S. Pat. No. 5,257,183, image reconstruction processing 18 begins by calculating derivative data from the acquired cone beam data, such data being calculated for those portions of integration planes Q which intersect both the scan path and the ROI. One such Q-plane is shown in FIG. 2. In the forenoted U.S. Pat. No. 5,257,183, Radon space driven conversion of the derivative data is used to develop an image reconstruction of object 16 on a display 20. However, image reconstruction processing 18 of the present invention uses a detector data driven technique, such as the simplified 3D backprojection technique described in the forenoted U.S. Pat. No. 5,881,123.

[0038] In accordance with the principles of the present invention, and simplified image reconstruction processing is provided, which processing does not require that the imaging apparatus provide circle scan paths near the top and bottoms of the ROI.

[0039] Furthermore, a radiation blocking shutter 22 is also provided in the FIG. 2 apparatus and selectively operated so as to reduce the radiation dose directed toward the object being imaged during that portion of the scan path which is proximate the top and bottom portions of the ROL More specifically, as will be described below, there are source positions proximate the top and bottom of the ROI where the cone beam data acquired on the top and bottom half of the detector, respectively, will not be used. Accordingly, at these source positions, shutter 22 is activated to block radiation directed to the portion of the detector that will not be used.

[0040] In order to more fully understand the principles of the present invention, some review of the previously noted U.S. patent application Ser. No. 09/343,770 will be provided. As was noted above, this technique for ROI imaging does not have the requirement that the imaging apparatus provide circle scans even with the top and bottom ends of the ROI. This is accomplished by processing the acquired cone beam data using a “local 2 ROI” approach, as was first described in the present inventor's U.S. patent application Ser. No. 09/052,415 entitled PRACTICAL CONE BEAM IMAGE RECONSTRUCTION USING LOCAL REGIONS-OF-INTEREST, filed Mar. 31, 1998, incorporated herein by reference, wherein specific sub-sets of Radon space data are developed. Each of the subsets of Radon space data is targeted for reconstructing a corresponding sub-section or “local” portion of a 2D ROI. Each local ROI is a 2D parallel projection of the ROI in the object on an individual one of the plurality of vertically oriented coaxial φ-planes that partition the Radon space. Due to the nature of a parallel projection, the upper and lower bounds of the local 2D ROI projection images do not suffer data corruption from x-rays emitted from the source which pass through areas outside the bounds of the local ROI. When a sufficient number of these local 2D ROI projection images are developed over 180 degrees (the sufficiency being determined by the desired resolution of the final image reconstruction), they are processed, using Radon inversion techniques, for directly developing contributions to the final image reconstruction, thereby obviating the requirement to provide circle scans at the top and bottom ends of the ROI.

[0041] In the forenoted U.S. patent application Ser. No. 09/343,770 these image reconstruction techniques are further simplified and extended for use in a Filtered Backprojection (FBP) image reconstruction arrangement.

[0042] The FBP image reconstruction technique was previously described in conjunction with FIGS. 1A and 1B.

[0043] When reconstructing an image using data acquired from source positions which are “internal” to the top and bottom 2π angular range of the scan path about the ROI, the simplified ramp filtering/data-combination masking technique known from the forenoted U.S. patent application Ser. No. 09/162,303 can be used. However, when reconstructing an image from data acquired at source positions which are near (i.e., “within”) the top and bottom 2π angular range of the scan path about the ROI, if there are no circle scans, in accordance with the teachings of the forenoted prior art U.S. patent application Ser. No. 09/343,770, the techniques shown in FIGS. 3A and 3B (and processed according to the forenoted FBP image reconstruction sub-steps (1)-(3) ) could be used. However, the 2D projections and 2D backprojections for each θ onto the detector plane for the views near the two ends of the spiral, is a very time-consuming task.

[0044] The present invention simplifies this task

[0045] More specifically, as shown in FIGS. 3A and 3B, in a set of cone beam data 300 acquired at a given source position S_i, a point C₀is marked on the x-axis of the detector. The technique for marling the point C₀is described in detail in the forenoted U.S. patent application Ser. No. 09/343,770, and will be described more briefly below in conjunction with FIG. 4. Next, an appropriate data combination mask 302 is applied to the cone beam data, and a plurality of lines L are formed in each set 300 of the acquired cone beam data in accordance with the forenoted sub-step (1) of FIG. 1A. As previously noted, the endpoints (and hence the length) of the line segments L are determined by the outer boundaries of the data combination mask for all source positions that are “internal” of the top and bottom 2π angular range in the scan path that are near the top and bottom edges, respectively, of the ROI. However, at source positions which are “within” the top and bottom 2π angular ranges (i.e., the outermost or first turns) in the scan path, the endpoints for at least some of the line segments drawn in the acquired cone beam data are determined by one boundary and the horizontal axis of the data combination mask, as shown in FIGS. 3A and 3B. That is, all line segments L which cross the horizontal axis to one side of a point C₀(which point is determined as shown in FIG. 4), have their endpoints determined by one of the upper or lower boundaries of the mask and its horizontal (i.e., x) axis, and all of those line segments crossing the horizontal axis to the other side of C₀, are either not used, or are determined by outer boundaries of the data combination mask. (Note, in the illustrated embodiment the x-axis of the data combination mask is coincident with the x-axis of the detector, but such coincidence is not required.)

[0046] The following analysis, in conjunction with FIG. 4, describes how to determine the point C₀in the data acquired at any source position S_i, scanning towards an end source position S₀. FIG. 4 illustrates a top view of a portion of the image reconstruction geometry, wherein a circle 402 is representative of the object and a circle 404 is representative of the spiral path traversed by the x-ray source S (shown as 10 in FIG. 2). The ending source position S₀is marked on scan path circle 404. As previously described, in order to provide an exact image reconstruction and avoid image artifacts, a parallel projection of the scan path must, at a minimum, begin and end outside of the local ROI of the object. Accordingly, a line 406 is drawn starting from source position S₀so as to be tangential to circle 404 (i.e., touching the outside edge of scan path 14). The Radon space φ-plane which forms an orthogonal intersection with line 406 defines the Radon space φ₀-plane, which is labeled 408.

[0047] Then, assuming an anti-clockwise rotation of the source S, mark a source position S_ion scan path circle 404. Next, one draws a line 410 which passes through source position S_iand is parallel with line 406, and therefore also bears an orthogonal relationship to the Radon space φ₀-plane 408. Then, C₀is defined as that point on the x-axis in the detector plane that acquires the cone beam data at source position S_i, which is intersected by line 410. (Note, the detector plane that acquired the cone beam data at source position S_i, i.e., detector plane 412, is determined by its orthogonal relationship to a line 414, which is drawn from source position S_ito the center of object circle 402.) This C₀is used, as noted above, to group the line segments L formed in the cone beam data during image reconstruction processing, in accordance with the manner of their length limitation. In this manner, the location of C₀for all of the source positions S can be determined. (The anti-clockwise rotation assumed in the above analysis can be easily modified for the opposite sense of source rotation.)

[0048] As described above, it is clear that the line integrals are not computed between fixed boundaries. Accordingly, the line integrals are not in the form of the projections of a function, and therefore the 2D step of the reconstruction processing (the forenoted FBP image reconstruction sub-steps (1)-(4)) cannot be simplified in the same way as was described in the forenoted U.S. Pat. No. 5,881,123.

[0049] The present invention simplifies the 2D step into 1D operations for the data acquired at the views near the two ends of the spiral.

[0050] In accordance with the principles of the present invention, these most time-consuming tasks are combined into convolving operations, and more specifically, 1D Hilbert Transform steps, thereby tremendously simplifying this portion of the image reconstruction processing.

[0051] In an effort to explain the theory behind the present invention, it is noted that the resulting 2D image 109 of FIG. 1B, after processing the cone beam data according to 4 FBP image reconstruction sub-steps described in the forenoted U.S. patent application Ser. No. 09/343,770, can be written as:

[0052] Y=D_t{∫dθ[B(θ)D_r(θ)P(θ)]}M(X)  (1)

[0053] where

[0054] X=cone beam image

[0055] M=masking

[0056] P(θ)=projection in the direction θ

[0057] D_r(θ)=differentiation w.r.t. the variable r for the projection at angle θ

[0058] B(θ)=backprojection in the direction θ

[0059] D_t=differentiation (spatial) in the projected scan path direction t

[0060] The projection and backprojection operations in Equation (1) are carried out for all θ in the angular range [θ_t−π/2, θ_t+π/2], where θ_tis the angular displacement (from the x axis) of the direction t, which direction is the projection of the current scan path tangent on the detector. The angular range [θ_t−π/2, θ_t+π/2] ensures that the product of the derivative d/dt with d/dr has the same sign for all the angles within the range; otherwise the signs of the product at different angles will be inconsistent.

[0061] For views in the “interior” of the spiral, the projection operation (1) is applied to the entire masked cone beam image M(X). It was shown in the forenoted U.S. Pat. No. 5,881,123 that under such circumstances, the 2D step (Equation (1)) can be greatly simplified, and the image Y (109 of FIG. 3B) can be written as:

[0062] Y=R_tM(X)+C  (2)

[0063] where R_tis the 1D ramp filtering operation in the direction t. The term C in Equation (2) is a 2D “correction” term obtained by convolving all the pixels intersecting the mask boundary with a space-variant 2D detector pixel point spread function representative of the effect of the image reconstruction processing upon that point.

[0064] For the views near the two ends of the spiral, however, as shown in FIGS. 3A and 3B, and described in the forenoted U.S. patent application Ser. No. 09/343,770, the projection operation (1) is applied to only part of the masked cone beam image M(X). Accordingly, the 2D reconstruction processing step cannot be simplified in the same way as was described in the forenoted U.S. Pat. No. 5,881,123.

[0065] As noted above, the present invention simplifies the 2D step (Equation (1)) into a series of 1D operations for the views near the two ends of the spiral.

[0066] As proof for the invention, it will first be shown that FBP image reconstruction sub-steps (1) through (3) of the above noted 4 FBP image reconstruction sub-steps can be combined into a 1D convolution step. As mentioned above, FBP image reconstruction sub-steps (1) through (3) are carried out for θ in the angular range [θ_t−π/2, θ_t+π/2], where θ_tis the angular displacement of the direction t which is the projection of the current scan path tangent on the detector. At each angle the direction of the variable r determines the direction of the derivative operation d/dr in sub-step (2) (thus the derivative operations in θ and θ+π will proceed in opposite directions). The angular range [θ_t−π/2, θ_t+π/2] ensures that the product of the derivative d/dt with d/dr has the same sign for all the angles within the range; otherwise the signs of the product at different angles will be inconsistent. To stress the importance of this angular range, we state it as a proposition:

[0067] Proposition 1.

[0068] The direction of the derivative operation d/dr in sub-step (2) follows the direction of the variable r in the angular range [θ_t−π/2, θ_t+π/2].

[0069] The 1D simplification of the 2D operations, as provided by the present invention, is based on the following analysis.

[0070] Theorem 1.

[0071] The combined operations of FBP image reconstruction sub-steps (1) through (3) in the angular range [θ_t−/2, θ_t+π/2] is equivalent to the 1D Hilbert transform in the direction t. Proof:

[0072] {∫dθ[B(θ)D_r(θ)P(θ)]}(f(ρ,α))=∫dθ∫dk(ik)F(k,θ)exp(ikρ cos(θ−α))=∫∫dθdk|k|H(k)F(k,θ)exp(ikρ cos(θ−α))  (3)

[0073] where r=ρcos(θ−α) H(k)=i sgn(k•t).

[0074] FIG. 5 spatially illustrates where the term “sgn(k•t)”, as used in equation (3), is +1 or −1.

[0075] Thus

[0076] {∫dθ[B(θ)(d/dr)P(θ)]}(f(ρ,α))=h*f

[0077] where $\begin{array}{cc}\begin{array}{c}h=\text{Fourier transform of }H\left(k\right)\\ =g\left(\nu \right)\delta \left(u\right)\end{array}& \left(4\right)\end{array}$

[0078] where (u, v), shown in FIG. 5, is the (x, y) coordinate system rotated such that the y′ axis coincides with the direction t. Accordingly:

[0079] u=x cos(θ_t−π/2)+y sin(θ_t−π/2)

[0080] v=−x sin(θ_tπ/2)+y cos(θ_t−π/2)

[0081] The 1D function g(v)=is the Fourier transform of the function i sgn(k). It's functional form is given by:

[0082] g(ν)=−1/(πν)

[0083] and is known as the Hilbert transform (well known to those skilled in CT image reconstruction technology, and described in detail, for example by R. N. Bracewell, in “The Fourier Transform and Its Applications”, 1978, McGraw-Hill, p.267).

[0084] Thus the result of FBP image reconstruction sub-steps (1) through (3) on the function f is the convolution of f with the point spread function (PSF) h. Even though h is a 2D function, Equation (4) indicates that it is zero everywhere except on the line θ=θ_t. Thus h is essentially a 1D function along the line θ=θ_t. The 1D function can be written as h(x,y;x′,y′;θ_t) to represent the value produced at a field point (x,y) as a result of applying the FBP image reconstruction sub-steps (1) through (3) on a source point (x′,y′) which has a density value of unity.

[0085] Theorem 2.

[0086] The combined operations of FBP image reconstruction sub-steps (1) through (3) in the angular range [θ₁, θ₂], where θ₁ε[θ_t−π/2, θ_t+π/2] and θ₂ε[θ_t−π/2, θ_t+π/2], is equivalent to the sum of two 1D Hilbert transforms: (½)h(x,y;x′,y′; θ₂−π/2)+(½)h(x,y;x′,y′; θ₁+π/2).

[0087] Proof of Theorem 2 is now provided, in conjunction with the illustrations in FIGS. 6A, 6B and 6C.

[0088] Employing a procedure similar to the one leading to Equation (3), we have

[0089] {∫_adθ[B(θ)D_r(θ)P(θ)]}(f(ρ,α))=H_A*f  (5)

[0090] where

[0091] h_A=Fourier transform of (i sgn_A(k))

[0092] where, as shown in FIG. 6A $\begin{array}{ccc}{\mathrm{sgn}}_A\left(k\right)& =1& \mathrm{k\epsilon A}\\ & =-1& \mathrm{k\epsilon A}+\pi \\ & =0& \text{Otherwise}\end{array}$

[0093] Considering two angular ranges A₁and A₂, shown in FIGS. 6B and 6C, respectively:

[0094] A₁=[θ₁,θ₁+π], t₁=direction vector at θ=θ₁+π/2, h₁=FT[i sgn(k•t₁)]

[0095] A₂=[θ₂−π,θ₂], t₂=direction vector at θ=θ₂−π/2, h₂=FT[i sgn(k•t₂)]

[0096] From inspection of FIGS. 6A-6C, is is clear that:

[0097] sgn_A(k)=(½)[sgn(k•t1)+sgn(k•t₂)] (i.e., FIG. 6A=½(FIG. 6B+FIG. 6C)

[0098] Thus

[0099] h_A=(½)[h₁+h₂)

[0100] Applying Threorem 1 to A₁and A₂yields the desired result.

[0101] Next, in accordance with the principles of the invention, Theorems 1 and 2 are applied to simplify the know prior image reconstruction processing for cone beam data which is acquired near the edges of the ROI, and without circle scans.

[0102] Without loss of generality, assume the spiral scanning motion is right-handed with the spiral axis pointing in the +y direction; the analysis can easily be modified for other scanning configurations. Thus the direction t, which is the projection of the current scan path tangent onto the detector, points in the direction θ_t[0,π/2].

[0103] We first consider the view angular range [π, 2π] measured from the spiral end, with the corresponding limits for the projection operation being as shown in FIG. 3B. Considering a source point P=(x′,y′), there are generally two possibilities:

[0104] Case 1.

[0105] Where y′<0 (i.e., all source points P are below the x-axis in the mask, i.e., are in its bottom half.)

[0106] Line integrals are computed on all line segments through P. The angular range of the projection and backprojection operations is therefore π. Application of Theorem 1 yields the result:

[0107] PSF at P=h(x,y;x′,y′;θ_t)

[0108] I.e., the result is a 1D Hilbert transformation in direction t for each point P on the detector which is below the x-axis.

[0109] FIG. 7 illustrates the taking of these 1D Hilbert transformations along a plurality of parallel line segments 700 drawn in the bottom half of the masked data.

[0110] Note, the results of these transformations, as well as the transformations noted below, are spread out over a virtual detector which is enlarged (compared with a detector having a height corresponding to the pitch of the scan path), so as to cover the entire ROI in the object. As previously noted, this enlargement is required because the calculated data affects the reconstruction of the entire Q plane, and not just the partial plane represented by the data combination mask.

[0111] Case 2. Where Y′>0

[0112] Let C₀be the point on the x axis which divides the x axis into data-combination limited and ROI-edge limited portions, as described in the forenoted U.S. Ser. No. 09/343,770. For the assumed scanning motion, it can be shown that the portion to the right of C₀is data-combination limited, and that the portion to the left of C₀is ROI-edge limited. Referring now to FIGS. 8A and 8B, illustrating a masked set of acquired cone beam data, let a point B=(x_b,y′) be a point on a horizontal line through point P such that X_b>x′, i.e. B is on the right hand side of P. Consider all the line integrals through the point P. These line integrals occur in the angular range bound by the line segment PC₀anti-clockwise through the line segment PB. Although not shown in FIGS. 8A and 8B, it is understood that the integration line segments are limited by the top and bottom boundaries of the data-combination mask as described in the forenoted U.S. Ser. No. 09/343,770. Let −θ₀(θ₀>0) be the angular displacement of the line segment PC₀from the x axis. Image reconstruction processing for line segments in this range of angular displacement resolves into one of two sub-cases.

[0113] Sub-Case 2.1, where θ_t−π/2>−θ₀+π/2 (i.e., a line drawn in the t direction crosses the x-axis to the right of C₀, as shown in FIG. 8A,) Consider the line integral on the line segment PB. Depending on the direction of integration, the line integral can be considered as belonging to either the projection at θ=π/2 or the projection at =−π/2. Since Proposition 1 requires that the projection/backprojection angle should be inside the angular range [θ_t−π/2, θ_t+π/2], θ=/2 is the projection/backprojection angle. Following this line of argument it can be shown that the projection/backprojection angular range corresponding to all the line integrals in the angular range bound by the line segment PC₀anti-clockwise through the line segment PB is A1∪A₂, where A₁=[θ_tπ/2, π/2]⊂[θ_t−π/2, θ_t+π/2] and A₂=[θ₀+3π/2, θ_t+π/2]. Application of Theorem 2 to A₁yields the result:

[0114] PSF at P (due to A₁)=(½)h(x,y;x′,y′;0)+(½)h(x,y;x′,y′;θ_t)  (6)

[0115] Application of Theorem 2 to A₂yields the result:

[0116] PSF at P (due to A₂)=(½)h(x,y′,x′,y′;θ₀)+(½)h(x,y;x′,y′;θ_t)  (7)

[0117] Combining the results (6) and (7) we have:

[0118] PSF at P=(½)h(x,y;x′,y′;+θ₀)+(½)h(x,y;x′,y′;0)+h(x,y;x′,y′;θ_t)

[0119] I.e., the result is the summation of three 1D Hilbert transformations for each point P which satisfies the criteria of this sub-case, as shown by lines 902, 904 and 906, respectively, in FIGS. 9A, 9B, and 9C.

[0120] Sub-Case 2.2. where θ_t−π/2 <−θ₀+π/2 (i.e.a line drawn in the t direction crosses the x-axis to the left of C₀. as shown in FIG. 8B) Following the same reasoning used above, it can be shown that the projection/backprojection angular range corresponding to all the line integrals in the angular range bound by the line segment PC₀anti-clockwise through the line segment PB is [π/2,θ₀+3π/2]⊂[θ_t−/2, θ_t+π/2]. Application of Theorem 2 yields the result:

[0121] PSF at P=(½)h(x,y;x′,y′;−θ₀+π)+(½)h(x,y;x′,y′; 0)

[0122] I.e., the result is the summation of two 1D Hilbert transformations for each point P which satisfies the criteria of this sub-case, each 1D transformation being in the direction as shown by lines 1002 and lines 1004, shown in FIGS. 10A and 10B, respectively.

[0123] Summarizing, for the view angular range [θ, 2π] measured from the end of the ROI, the above analysis has shown that applying FBP image reconstruction sub-steps (1) through (3) to any point on the detector, results in a number of essentially 1D PSF h, which, in accordance with the present invention, are additively combined for constructing the final backprojection image Y (109 of FIG. 1B), as follows:

[0124] For points in that portion of the image which is below the x axis (case 1):

[0125] 1D Hilbert transform in the direction of the scan path [h(x,y;x′,y′;θ_t)] (as shown by FIG. 7).

[0126] For points in that portion of the image which is above the x axis (case 2):

[0127] 1D Hilbert transform in the x direction [h(x,y;x′,y′;0)] (as shown by the combination of FIGS. 9B and 10B).

[0128] Determine the point C₀on the x axis where the current source position projects in the direction which projects the end source position onto the edge of the spiral cylinder.

[0129] 1D Hilbert transform in the direction of the scan path the portion of the image above the x axis and to the right of the line (infinitely extended) through C₀in the direction of the scan path fh(x,y;x′,y′;θ_t)] (as shown by FIG. 9C).

[0130] 1D Hilbert transform along each line intersecting the point C₀[h (x,y;x′,y′;−θ₀) or h(x,y;x′,y′;θ₀+π)], respectively, as shown by the FIGS. 9A and 10A).

[0131] In a similar fashion one can work out the PSF for cone beam images acquired during the 0 to π angular range of the scan path about the top end of the ROI, as follows:

[0132] Case 1. y′>0

[0133] PSF at P=0

[0134] Case 2. y′<0

[0135] Case 2.1 θ_t−π/2 >−θ₀+π/2 (Similar to FIG. 8A but with the location of P and C₀interchanged)

[0136] PSF at P=(½)h(x,y;x′,y′;−θ₀+π)+(½)h(x,y;x′,y′;π)

[0137] Case 2.2 θ_t−π/2<−θ₀+π/2 (Similar to FIG. 8B but with the location of P and C₀interchanged)

[0138] PSF at P=(½)h(x,y;x′,y′;−θ₀)+(½)h(x,y;x′,y′;π)+h(x,y;x′,y′;θ_t)

[0139] Summarizing, for the view angular range [0,π] measured from the spiral end we have shown that applying sub-steps (1) through (3) to any point on the detector results in a number of essentially 1D PSF h:

[0140] For the portion of the image above the x axis (Case 1. y′>0):

[0141] No computation 16 For the portion of the image below the x axis (Case 2. y′<0):

[0142] 1 D Hilbert transform in the −x direction [h(x,y;x′,y′;π)].

[0143] Determine the point C₀on the x axis where the current source position projects in the direction which projects the start source position onto the edge of the spiral cylinder.

[0144] 1 D Hilbert transform in the direction of the scan path the portion of the image below the x axis and to the right of the line (infinitely extended) L_cthrough C₀in the direction of the scan path [h(x,y;x′,y′;θ_t)].

[0145] 1 D Hilbert transform along each line intersecting the point C₀[h(x,y;x′,y′;−θ₀) and h(x,y;x′,y′;−θ₀+π)] respectively for the portions to the right and the left of the line L_c.

[0146] Since the above processing corresponds to the FBP substeps (1)-(3), it is still necessary to do sub-step (4), to complete the 2D processing steps. Since sub-step (4) is also a 1D operation, i.e. differentiating the backprojection image in the projected scan path direction t, the entire 2D step represented by Equation (1), in accordance with the present invention, is reduced to a number of 1D operations for the views near the 2 ends of the spiral scan, i.e. views in the first and the last 2π angular range.

[0147] For the cone beam images acquired near the bottom of the ROI, i.e., the start of the scan path, the image processing is substantially the same as described above for processing cone beam images acquired near the top of the ROI, except that the sets of line segments shown FIGS. 3A and 3B that are bound by the x-axis and the bottom boundary of the mask, are now bound by the x-axis and the top boundary of the mask. This analysis can proceed in a similar manner to make corresponding modifications to the above description.

[0148] Thus, there has been shown and described a novel 3D CT imaging method and apparatus for obtaining an exact image reconstruction which eliminates the need for performing 2D projections and 2D backprojections for each θ onto the detector plane for the views of the ROI near the two ends of the spiral scan path, thereby greatly simplifying this portion of the FBP image reconstruction processing. Many changes, modifications, variations and other uses and applications of the subject invention will, however, become apparent to those skilled in the art after considering this specification and its accompanying drawings, which disclose a preferred embodiment thereof. For example, many alternative arrangements were noted in the above text. For example, the pitch of the spiral scan could be variable, instead of fixed, and the convolution processing could be done in the spatial frequency domain using the Fourier transform of the Hilbert kernal. Furthermore, although in the preferred embodiment the x-axis of the detector was coincident with the x-axis of the data-combination masks, such coincidence is not required, nor is it required to provide a complete blocking of the radiation exposure to one side of the axis (even a partial blocking of the radiation to one side of the axis would be desirable). All such changes, modifications, variations and other uses and applications which do not depart from the teachings herein are deemed to be covered by this patent, which is limited only by the claims which follow as interpreted in light of the foregoing description.