The following relates to the diagnostic imaging arts. It finds particular application in helical conebeam computed tomography imaging, and will be described with particular reference thereto. However, it also finds application in other types of tomographic imaging.
Single-slice x-ray transmission computed tomography imaging systems employ limited beam dimensions, such as an essentially planar fan-beam that is rotated in a circular orbit around an imaging subject. The acquired projection data is reconstructed to generate an image slice. In multi-slice fan-beam imaging, data are acquired for a plurality of axially spaced-apart circular orbits, and reconstructed to produce a corresponding plurality of axially spaced-apart image slices. Multi-slice computed tomography imaging is conducive to fast and accurate image reconstruction due to the simplified fan-beam geometry. In one suitable reconstruction method, each slice is reconstructed separately using convolution backprojection. In another reconstruction method, the fan-beam projection data are rebinned into substantially parallel-ray views that are backprojected to reconstruct the image. Even in multi-slice methods, the rate of data acquisition is limited due to the small volume sampled by the essentially planar fan-beam and the discontinuous axially stepped orbiting of the x-ray source.
Helical conebeam computed tomography was developed to overcome these data acquisition rate constraints. A continuous helical x-ray source orbit is combined with a conebeam geometry that interacts with a three-dimensional volume to enable much higher data acquisition rates. A two-dimensional array of x-ray detector elements is arranged to receive the conebeam after it passes through an examination region in which the imaging subject is disposed. The detector array can be flat (flat detector geometry), curved (curved detector geometry), or otherwise shaped.
A disadvantage of helical conebeam computed tomography is that reconstruction of the conebeam projection data into an image representation is computationally complex. A given voxel of the image volume is sampled by non-coplanar projections spanning a portion of the helical trajectory of the x-ray source. In contrast, with the fan-beam geometry projections for a given voxel are substantially coplanar, that is, lie in a single plane.
Various conebeam reconstruction methods have been developed. Some of these techniques, such as the N-PI method, are not exact, in that they fail to accurately account for the cone angle. These inexact reconstruction methods introduce artifacts into the reconstructed image because they do not accurately account for the cone angles of the projections, that is, angular deviations of the projections from the axial plane. By employing various approximations relating to the cone angle, inexact reconstruction methods can provide rapid reconstruction with acceptably limited image artifacts for data acquired using small cone angles. However, these inexact reconstructed imnages become increasingly degraded by artifacts as the cone angle increases. Conebeams with large cone angles provide large beam interaction volumes and advantageously enable fast data acquisition.
For data acquired with large cone angles, an exact reconstruction method is preferably used that accurately accounts for projection components in the cone angle direction.
Prior exact reconstruction methods are more computationally intense than inexact methods due to the more complex conebeam geometry. An exact conebeam reconstruction method has been developed by Katsevich (see for example Katsevich et al, Proceedings SPIE Medical Imaging Conference, San Diego, Calif. (February 2003)).
However, this exact conebeam reconstruction method employs a voxel-based coordinate system that does not comport with conebeam data acquisition detector geometries.
Conversion of acquired projection data, which are defined using coordinates such as x-ray detector row and column values, projection angle coordinates, helix angle, or the like, into a detector-independent voxel-based coordinate system is computationally challenging and involves additional rebinning operations that substantially increase reconstruction time. Moreover, the previous exact conebeam reconstruction method is not readily adapted to existing backprojector processing pipelines such as are commonly employed in reconstructing parallel-rebinned projection data.
The present invention contemplates an improved apparatus and method that overcomes the aforementioned limitations and others.
According to one aspect, a conebeam computed tomography imaging system is disclosed. An x-ray source produces an x-ray conebeam directed into an examination region. The x-ray source is arranged to traverse a generally helical trajectory around the examination region.
An x-ray detector array is arranged to detect the x-ray conebeam after passing through the examination region. The x-ray detector array generates projection data in native scan coordinates defined with reference to the detector array. An exact reconstruction processor performs an exact reconstructing of conebeam projection data produced by the detector array into an image representation. The reconstructing is performed in the native scan coordinates
According to another aspect, a conebeam computed tomography imaging system is disclosed that produces conebeam computed tomography projection data having native scan coordinates that include at least a helix angle, a projection fan coordinate, and a projection cone coordinate. A means is provided for computing filtered projection data. The means for computing includes a differentiating means for computing a derivative of projection data with respect to the helix angle at fixed projection direction and a convolving means for convolving projection data with a kernel function. The convolving is performed in the native scan coordinates. A means is provided for backprojecting the filtered projection data to obtain an image representation.
According to another yet aspect, an exact reconstruction method is provided for reconstructing conebeam computed tomography projection data having native scan coordinates that include at least a helix angle, a projection fan coordinate, and a projection cone coordinate. Filtered projection data is computed by a combination of computing a derivative of projection data with respect to the helix angle at fixed projection direction and convolving projection data with a kernel function. The convolving is performed in the native scan coordinates. The filtered projection data is backprojected to obtain an image representation.
One advantage resides in fast, exact image reconstruction of conebeam imaging data acquired using a large cone angle.
Another advantage resides in simplified image reconstruction computations by obviating rebitiing operations associated with transforming tomographic projection data into a detector-independent voxel-based coordinate system.
Another advantage resides in optional incorporation of a backprojector pipeline configured for processing parallel-rebinned projection data in performing exact reconstruction of conebeam projection data.
Yet another advantage resides in straightforward accommodation of various detector geometries such as flat detector geometries and curved detector x-ray source-focused detector geometries.
Numerous additional advantages and benefits will become apparent to those of ordinary skill in the art upon reading the following detailed description of the preferred embodiments.
The invention may take form in various components and arrangements of components, and in various process operations and arrangements of process operations. The drawings are only for the purpose of illustrating preferred embodiments and are not to be construed as limiting the invention.
With reference to
To effect a helical trajectory of the x-ray source 12 about an imaging subject, the imaging subject is placed on a couch 20 or other support. The couch moves linearly along a z-direction as indicated. The x-ray source 12 and the x-ray detector 16 are oppositely mounted respective to the examination region 14 on a rotating gantry 22, such that rotation of the gantry 22 effects rotation of the x-ray source 12. Rotation of the gantry 22 along with simultaneous, continuous linear motion of the couch 20 effects a helical trajectory of the x-ray source 12 around the imaging subject disposed on the couch 20.
The x-ray detector 16 is shown mounted on the rotating gantry 22 such that it rotates along with the x-ray source 12 to intercept the x-ray conebeam throughout the helical trajectory. However, it is also contemplated to replace the x-ray detector 16 by an x-ray detector band mounted around a stationary gantry 24.
In operation, during helical orbiting of the x-ray source 12 relative to the imaging subject, the x-ray conebeam is projected into the examination region 14 where it interacts with the imaging subject. Some portion of the x-rays are absorbed by the imaging subject to produce a generally spatially varying attenuation of the x-ray conebeam. The x-ray detector 16 measures the x-ray intensities across the conebeam to generate x-ray absorption data that is stored in a projection data memory 30.
With continuing reference to
The (x,y,z) stationary coordinate system is a conventional Cartesian coordinate system and is provided for reference. The helix angle X identifies an angular position of the x-ray source 12 along the helical orbit.
With particular reference to
A vertex position A(λ) of the x-ray source 12 along the helical trajectory is given by helix angle λ for the curved detector geometry. In the stationary (x, y, z) coordinate system, the vertex points a(λ) of the helix, that is, the trajectory of the x-ray source 12, is given by:
where Ro is a source distance, that is, a distance from the x-ray source 12 to the helix axis, P is the helical pitch, that is, the axial or z-distance the x-ray source moves relative to the imaging subject over a single helical trajectory turn, and λo and zo are a reference helix angle and reference axial position, respectively.
With particular reference to
The coordinate w is parallel. to the axial or z-direction. The source-to-detector distance D of the curved geometry is also defined for the flat detector geometry. Equation (1) specifying the vertex position a(λ) of the x-ray source 12 in the stationary (x,y,z) coordinate system also applies to the flat detector geometry of
The detector arrays 16c, 16f are shown for convenience as having a small number of detector elements. In
With returning reference to
Optionally, the filtered backprojector processor 42 includes weighting of the projection data with respect to detector position (aperture weighting) to reduce artifacts introduced by large projection components in the cone angle direction. The filtered backprojector processor 42 can also include weighting of the projection data with respect to the helix angle (angular weighting) to smooth the data at angular discontinuities, to combine angularly complementary projection data, or the like. Suitable three-dimensional filtered backprojector pipes are described, for example, in U.S. Pat. No. 6,104,775 issued to Tuy, and in U.S. patent application Ser. No. 10/274,816 by Heuscher et al. filed on Oct. 21, 2002.
The resultant image representation is suitably processed by a video processor 50 to generate a three-dimensional rendering, one or more image slices, or other visual representation of the reconstructed image that is displayed on a user interface 52. Rather than a video display, the image representation can be formatted by a printer driver and printed out using a printer, transmitted over an electronic network, stored electronically, or otherwise utilized. Preferably, the user interface 52 communicates with a computed tomography controller 54 to enable a radiologist or other operator to initiate imaging or otherwise control operation of the computed tomography scanner 10.
The inexact image reconstruction processing employing the parallel rebinning processor 34, the ramp convolution processor 36, and the parallel three-dimensional backprojector processor 42 may be adequate for certain applications. However, if an exact reconstruction which fulfills all the requirements of the three-dimensional Radon transform is desired, then an exact filtered backprojection conebeam reconstruction process with a one-dimensional shift-invariant filtering is employed. The disclosed exact reconstruction method is compatible with projection data that is axially or z-truncated. When performing exact conebeam reconstruction, a different processing path starting at a finite derivative processor 60 is engaged. That is, rather than inputting the projection data into the parallel rebinning processor 34, the data is input into the finite derivative processor 60.
With continuing reference to
A post-cosine weighting processor 80 performs a cos(α) weighting in the case of the curved detector geometry of
The image representation is suitably processed by the video processor 50 and displayed on the user interface 52, or otherwise utilized.
Alternatively, the data output by the convolution processor 64 is input to an inverse cosine weighting processor 80′, parallel-rebinned by a parallel rebinning processor 34′, and backprojected by the 3D parallel backprojection processor 42 to generate an image representation that is stored in the image memory 44. For the flat detector geometry of
Those skilled in the art will appreciate that the described hybrid convolution processing operates entirely within the native scan coordinates of the curved or flat detector 16c, 16f. This increases reconstruction speed versus exact reconstruction performed in a detector independent coordinate system by eliminating computationally intensive rebinning operations associated with converting to the detector-independent voxel-based geometry. Moreover, reconstruction accuracy is promoted by avoiding interpolations involved in the rebinning to a voxel-based coordinate system.
The reconstruction processing components can be physically implemented in various ways. Some or all components can be software modules that are executed on one or more computers. Some or all components can be application-specific hardware pipeline components. Some or all of the reconstruction processing components can be integrated into the user interface 52, and/or into the computed tomography scanner 10, or can be stand-alone components. Hardware pipeline components are readily embodied as computer cards that insert into interface slots of a computer to enable integration of hardware- and software-based reconstruction processing components. Moreover, a given component, such as the backprojector 42, 82 can be implemented partially in software and partially as a hardware pipeline. Those skilled in the art can readily construct other physical implementations using other combinations of hardware and software.
With reference to
A κ-plane is a plane that has three intersections with the helix such that the middle helix intersection point is equally angularly spaced in helix angle λ from the two other helix intersection points. There exists a family of K,-planes corresponding to each vertex point a(λ). These corresponding κ-planes are indexed by an angle W that lies in a range (−π, π). The κ-plane of angle ψ at the vertex point a(λ) is designated herein as K(λ,ψ) and contains the vertex points a(λ), a(λ+ψ), and a(λ+2ψ). A unit vector normal to K(λ,ψ) and having an acute angle with respect to the helix axis is designated herein as n(λ,ψ), and is given by:
where the symbol “x” denotes a cross-product. As the angle ψ tends toward zero, the κ-plane K(λk,ψ) converges to a plane that is tangent to the helix H at a(λ) and is parallel to the detector v coordinate, that is:
All κ-planes have the following property, described with reference to
With returning reference to
eu(λ)=[−sin(λ+λo), cos(λ+λo), 0]ev(λ)=[−cos(λ+λo), −sin(λ+λo), 0]ew=[0, 0, 1]
where the directions u, v, w are set forth in Table I. In this rotated coordinate system and for the curved detector geometry of
Conversely, given the projection direction vector θc the curved detector coordinates (α,w) can be computed as:
where the dot operator denotes a dot product. Using these coordinates, the derivative of g(λ,α,w) with respect to helix angle λ at fixed projection direction θc for the curved detector coordinate system is given by:
The derivative of Equation (7) is preferably implemented by the finite derivative processor 60 using a discrete finite difference derivative computation that arithmetically combines four neighboring projection samples. For an exemplary case where the sampling increment in coordinate α is one-fourth the sampling increment in the helix angle λ, that is, Δλ=4Δα, a suitable finite derivative is:
D′nj=a0Dn−1j+a1Dnj+a2Dnj+1+a3Dn+1j+1
where n is the discrete projection sample index along the helix angle direction λ, j is the discrete projection sample index along the a coordinate direction, and the constants a0, a1, a2, a3 are given by:
a0=⅛
a1=⅞
a2=−⅞
a3=−⅛
A half-sample shift along the helix angle λ is introduced by the finite derivative of Equations (8) and (9).
It will be recognized that the finite derivative of Equation (8) is advantageously performed on the projection data in the curved detector coordinate system without computationally intensive rebinning. Equations (8) and (9) are exemplary equations for the specific case of finite derivative computation using four neighboring projections with a spacing ratio of Δλ=4Δα. Those skilled in the art can readily construct other finite derivatives using another number of neighboring projections, for other spacing ratios, and so forth.
With reference to
The differentiated and length-normalized projection data output by the cone angle length correction processor 62 is processed by the convolution processor 64. Preferably, the convolution is implemented by the FFT convolution processor 72. The convolution is a one-dimensional convolution with respect to coordinate a performed along intersections of κ-planes K(λ,ψ) with the curved detector 16c. In other words, the convolution is a one-dimensional convolution with respect to coordinate a performed at a fixed-angle ψ. Hence, prior to input to the FFT convolution processor 72, the differentiated and length-normalized projection data is rebinned by the forward height rebinning processor 70 to get constant ψ surfaces according to:
g3(λ,α,ψ)=g2(λ,α,wκ(α,ψ)) (11)
over all ψ in a range [−π/2-αm, π/2+αm] where αm is a fan angle defined by the size R of the field of view and the helix radius Ro, that is, αm=arcsin(R/Ro), with wκ(α,ψ) given by:
At a fixed angle ψ, it will be recognized that Equation (12) describes a curve in the detector area that is the intersection between the detector array of the curved detector 16c and the κ-plane K(λ,ψ). This curve is referred to herein as a κ-curve of angle ψ.
With reference to
The FFT convolution processor 72 is applied to the rebinned data g3(λ,α,ψ) as a one-dimensional convolution in angle α for constant angle ψ according to:
where hH is a kernel based on a Hilbert transform. A suitable kernel hH is abs(sin(x))/sin(x) which in the s-domain is:
However, in view of the half-sample shift introduced by the selected finite derivative of Equations (8) and (9), a compensating half-sample shift is preferably incorporated into the kernel hH.
After the one-dimensional convolution processor 72, the projection data g4(λ,α,ψ) is processed by the backward height rebinning processor 74 to obtain rebinned projection data g5(λ,α,w) according to:
g5(λ,α,w)=g4(λ,α,ψmin(α,w)) (15)
where ψmin is the angle ψ of smallest absolute value that satisfies the equation:
Reviewing the operation of the exemplary convolution processor 64, it will be seen that the output g5(λ,α,w) at a given detector location (αo,w0) is obtained by convolving projection data g2(λ,α,w) on the κ-curve of smallest absolute ψ value that passes through (αo,wo).
For processing data in the source-focused curved detector geometry, the output of the convolution processor 64 is preferably processed by the post-cosine weighting processor 80. This weighting adjusts for the convolution in the fan direction according to:
gF(c)(λ,α,w)=cos(α)g5(λ,α,w) (17)
to produce the final filtered data gF(c)(λ,α,w).
The filtered data gF(c) is suitable for subsequent processing by a conebeam backprojector processor 82 to produce an exact reconstruction of conebeam projection data. In a suitable embodiment, the conebeam backprojector processor 82 computes the image f(x) according to:
where λi(x) and λo(x) are the endpoints of the π-line through x with λi(x)<λo(x), and:
The backprojection of Equations (18)-(21) implemented by the conebeam backprojector processor 82 is exemplary only. Those skilled in the art can employ other suitable backprojection processes. For instance, in an alternative backprojection path shown in
In performing above-described exact reconstruction method, the number of detector rows should be large enough so that the filtered data can be computed for points of the region of interest. The number of detector rows, the helix pitch P, the helix radius Ro, and the size of the field of view can be interrelated. One example of a suitable detector region for conebeam reconstruction by the backprojection of Equation (18) is known in the art as the Tam-Danielsson window. An exemplary Tam-Danielsson window is shown as the shaded area on the detector in
and a maximum value wwp(a) given by:
The convolution processor 64 is a non-local process that uses projection data outside the Tam-Danielsson window. However, as shown in
where dw is a thickness of the detector rows.
The exact conebeam reconstruction method described with reference to Equation (4) through Equation (24) is suitable for reconstructing projection data formatted in the source-centered curved detector geometry of
With reference to
The projection direction vector θf for the flat detector geometry is given by:
Conversely, given the projection direction vector θf the flat detector coordinates (u,w) can be computed as:
where the dot operator denotes a dot product. Using these coordinates, the derivative of projection data g(λ,u,w) for the flat detector coordinate system with respect to helix angle λ at fixed projection direction θf is computed by the finite derivative processor 60 in accordance with Equations (7)-(9), making the substitutions of θf for θc and u for α in Equation (7), to obtain differentiated projection data g1(λ,u,w)=g′(λ,θf).
The length correction in the case of the flat detector geometry is computed by the cone angle length correction processor 62 as:
which is analogous to the normalization of Equation (10) for the curved detector geometry.
The forward height rebinng processor 70 performs rebinning of the flat detector geometry projection data according to:
g3(λ,u,ψ)=g2(λ,u,wκ(u,ψ)) (28),
where wκ(u,ψ) is:
Equations (28) and (29) are analogous to Equations (11) and (12) of the curved detector geometry processing. Analogously to the case of the curved detector geometry, Equation (29) defines κ-curves K(λ,ψ) in the flat detector geometry that correspond to intersections of κ-planes with the flat detector 16f.
The one-dimensional convolution processor 72 performs the following one-dimensional convolution with respect to coordinate u at constant angle ψ:
where a suitable kernel hH is given in Equation (14). The reverse height rebinning processor 74 rebins the convolved projection data g4(λ,u,ψ) to produce filtered projection data gF(f)(λ,u,w) according to:
gF(λ,u,w)=g4(λ,u,ψmin(u,w)) (31),
where ψmin is the angle of smallest absolute value that satisfies:
For the flat detector geometry, the post cosine weighting processor 80 or the inverse cosine weighting processor 80′ is bypassed, and so the output of the reverse height rebinning processor 74 is the filtered projection data gF(f)(λ,u,w). It will be appreciated that the processes set forth in Equations (28)-(32) obtain gF(f)(λ,uo,wo) from the projection data g2(λ,u,w) by convolving along the κ-curve of smallest absolute angle ψ that passes through (uo,wo).
The filtered data gF(f) is suitable for subsequent processing by the conebeam backprojector processor 82 to produce an exact reconstruction of conebeam projection data. In a suitable embodiment, the conebeam backprojector processor 82 computes the image f(x) according to:
analogous to the backprojecting of the filtered projection data gF(c) of the curved coordinate geometry given in Equation (18). The parameter v*(λ,x) is the same as that given in Equation (19) for the curved detector geometry, while the parameters u*(λ,x) and w*(λ,x) are given in the flat detector geometry by:
As in the case of backprojecting the filtered projection data gF(c) of the curved coordinate geometry, the backprojection process of Equations (19) and (33)-(35) is exemplary only. In one alternative backprojection approach, the filtered backprojection data gF(f) is input into the parallel rebinning processor 34′, and the parallel-rebinned filtered backprojection data is backprojected using the 3D parallel backprojector 42.
In the case of the flat detector geometry, a suitable number of detector rows Nrows is estimated as follows. The Tam-Danielsson window is the set of (u,w) points with u in a range [−um,um] where um=D tan(αm)=D tan(arcsin(R/Ro)), and w in a range [wbottom(u),wtop(u)], where:
Computation of the filtered data at any point in the Tam-Danielsson window employs data on κ-curves that intersect the Tam-Danielsson window. This data generally corresponds to data acquired using detector rows between ubottom and utop inclusive, and so a suitable number of detector rows is about:
which can be rewritten as:
where dw is the thiclkess of the detector rows.
The exact nature of the exact reconstructions described herein has been demonstrated by showing equivalency of image representations reconstructed using the techniques disclosed herein and image representations reconstructed using the exact method of Katsevich. The exact reconstruction methods described herein advantageously are performed entirely within the flat or curved detector coordinate system. Computationally intensive transformations of projection data into a detector-independent coordinate system is thus obviated, providing a more rapid reconstruction. The exact reconstruction methods described herein optionally incorporate the three-dimensional backprojector 42. This enables the exact reconstruction to take advantage of three-dimensional backprojector capabilities such as complementary combination of redundant projection data..
The invention has been described with reference to the preferred embodinents. Obviously, modifications and alterations will occur to others upon reading and understanding the preceding detailed description. It is intended that the invention be construed as including all such modifications and alterations insofar as they come within the scope of the appended claims or the equivalents thereof.
Filing Document | Filing Date | Country | Kind | 371c Date |
---|---|---|---|---|
PCT/IB04/00386 | 2/9/2004 | WO | 8/10/2005 |
Number | Date | Country | |
---|---|---|---|
60447426 | Feb 2003 | US | |
60482380 | Jun 2003 | US |