1. Field of the Invention
The present invention concerns a method for reconstruction of computed tomography (CT) representations from x-ray CT data sets of an examination subject, of the type wherein scanning of the examination subject ensues on a spiral helical path and, for reconstruction, a differential backprojection is implemented, followed by a Hilbert transformation over a surface formed by π-lines.
2. Description of the Prior Art
Differential backprojection (DBP) followed by a Hilbert transformation is generally known in connection with a CT spiral scan for reconstruction of computed tomography representations from x-ray CT data sets of an examination subject along lines known as π-lines, which are straight lines that intersect the spiral path twice at a distance of less than one full revolution. Such reconstruction is described in J. Pack, F. Noo and R. Clackdoyle, “Cone-beam reconstruction using the backprojection of locally filtered projections”, IEEE Trans. Med. Imag., vol. 24, no. 1, pp. 70-85, January 2005. In the known methods, the reconstruction always occurs along π-lines that lie on actual measured beam positions. This makes the method relatively inflexible.
An object of the present invention is to provide an improved and more flexible method for reconstruction of CT representations making use of π-lines.
The method, CT apparatus and computer-readable medium in accordance with the invention are based on the insight that it is advantageous when the conventional limitation of the conducted reconstruction to actually measured beam paths is abolished. According to the invention, this can occur by arbitrary virtual beams being interpolated with the use of actual measured, adjacent beams.
Therefore, in accordance with the invention, in a method for reconstruction of computer tomography representations from x-ray CT data sets of an examination subject, scanning of the examination subject ensues on a spiral path and, for reconstruction, a differential backprojection (DBP) is conducted followed by a Hilbert transformation over a surface formed by π-lines, and interpolated detector data are used for the reconstruction.
The spiral path is provided by a(λ)=[R0 cos(λ+λ0),R0 sin(λ+λ0),z0+hλ], wherein λ serves as a free index, R0 stands for the spiral path radius and h represents the table feed per spiral revolution. λ0 and z0 indicate the position of the x-ray source for the index λ=0.
It is advantageous when the detector data are interpolated between actual measured detector data, such that the associated π-lines are combined into a number of surfaces on which the π-lines respectively appear in parallel, projected on a plane perpendicular to a z-axis (=system axis) of a CT system used for examination. The interpolated detector data can be selected to form the π-lines, such that the projections of the π-lines are equidistantly formed on the (x,y) plane in a Cartesian (x,y,z) coordinate system. It should be noted that π-lines whose projections lie parallel to the (x,y) plane do not, by design, run parallel relative to the z-axis.
The individual surfaces that are formed from π-lines in this manner are indexed across the central π-line on this surface (thus the π-line that intersects the z-axis), and in fact across the first intersection point a(λfilt) of the source trajectory with this π-line. Each λfilt thus uniquely designates a surface of π-lines.
To reconstruct a volume described in classical (x,y,z) coordinates, it is necessary to initially cover this volume with surfaces of therefore π-lines over a sufficient range of filt.
Two fundamental preferred embodiments are possible in the method according to the invention for reconstruction of the object function on these surfaces, namely:
(A) Backprojection is conducted on a (s,τ)-grid rotated with λfilt in a cylindrical (s,τ,λfilt)-coordinate grid, wherein the final reconstruction occurs via application of the inverse Hilbert transformation in the same geometry, and the results are subsequently interpolated on a Cartesian (x,y,z) coordinate grid with:
s=−x sin(λfilt+λ0)+y cos(λfilt+λ0)
τ=−x cos(λfilt+λ0)−y sin(λfilt+λ0)
x,y,z=the Cartesian coordinates,
wherein the following relationship exists between the coordinates:
and
R0 is the radius of the spiral path
or
(B) the backprojection is executed on the surface of the π-lines across an (x,y)-grid in an (x,y, λfilt)-coordinate grid, the backprojection results are interpolated in a cylindrical (s,τ, λfilt)-coordinate grid in order to implement the inverse Hilbert transformation, and the results are subsequently interpolated on a Cartesian (x,y,z) coordinate grid with:
s=−x sin(λfilt+λ0)+y cos(λfilt+λ0)
τ=−x cos(λfilt+λ0)−y sin(λfilt+λ0)
x,y,z=the Cartesian coordinates,
wherein the following relationship likewise exists between the coordinates:
and
R0 is the radius of the spiral path.
In both embodiments A) and B), three different variants can be used. These are:
(a) a variant in which the derivations necessary for the DBP are effected exclusively in detector coordinates,
(b) a variant in which the DBP is implemented with the aid of a derivation according to the source position A given a fixed beam direction in space,
(c) a variant in which a rebinning in the pseudo-parallel rebinning geometry (=“wedge” geometry) that is defined by the following equations
initially occurs before the backprojection. The derivation and backprojection can subsequently be conducted in this geometry.
In the following the invention is described in detail with the use of
An implementation procedure is described for an effective two-step Hilbert reconstruction employing a differential backprojection (DBP) followed by an inverse Hilbert transformation (HT) on π-lines. This procedure is adapted to data sets that originate from curved detectors and are based on a reconstruction along theoretical π-lines. By theoretical π-lines, what is understood is a subset of the infinite set of π-lines which would be measured in a continuous measurement system. In contrast to real π-lines which exclusively represent those lines in space along which actual projections were obtained. This selection enables a free variation of the reconstruction grid within the x,y-plane.
Six different variations are within the scope of this implementation procedure. In the following, these six variants of the algorithm are described with their advantages and disadvantages.
The revolution path is provided by the equation
a(λ)=[R0 cos(λ+λ0),R0 sin(λ+λ0),z0+hλ], EQ 1
wherein λ describes the rotation angle of the source in the interval [0,λmax], R0 describes the spiral radius and 2πh describes the spiral advance. The revolution path is determined by λ0 and z0 such that, at λ=0, the source is positioned at the angle λ0 in the plane z=z0.
In contrast to the standard (x,y,z)-geometry, a rotating coordinate system is thus used which moves with the detector system corresponding to
e
u(λ)=[−sin(λ+λ0),cos(λ+λ0),0] EQ. 2
ev=[0,0,1] EQ. 3
e
w(λ)=[cos(λ+λ0),sin(λ+λ0),0]. EQ. 4
The detector consists of a field of Nrows×Ncols elements which are arranged in columns parallel to the unit vector ev and form arcs line-by-line through the spiral point around a line La parallel to the z-axis. It is thereby counted in the direction of the unit vectors eu and ew, wherein the detector units are described by the angle γ and the detector row index w while the detector does not completely lie in the plane that is formed by the unit vectors eu and ev. Furthermore, the distance from the line La through the focus towards the detector is designated with D, and the density of the three-dimensional examination subject is designated with f(x).
In the present notation, the measurements are described with:
A change of the standard (x,y,z)-geometry is required for the backprojection so that sets of parallel, theoretical π-lines can be worked with. We limit ourselves (without limiting the generality) to π-lines which exhibit a positive slope and form surfaces of π-lines that cover the desired reconstruction volume.
A surface of π-lines (as is shown in
For indexing on each π-line surface, s,τ-coordinates are used that are obtained by rotating the x- and y-coordinates corresponding to Equations 8 and 9:
e
s(λfilt)=[−sin(λfilt+λ0),cos(λfilt+λ0),0] EQ. 8
e
r(λfilt)=[−cos(λfilt+λ0),sin(λfilt+λ0),0]. EQ. 9
In other words, s describes the distance of projections of the π-lines on the x-y-plane from the origin and τ is a coordinate along these projections. Furthermore, a variable t along the π-lines is introduced, wherein the projection of t on the x-y-plane corresponds to the variable τ. The z-position at a point indexed by s,τ,λfilt is thus provided by:
For a volume V to be reconstructed in Cartesian coordinates, this equation allows us to determine the range of λfilt over which the backprojection should be executed in order to completely cover the volume V. The distance between the surfaces of the π-lines over which the backprojection is executed is then described with:
Δλfilt=Δz/h, EQ. 11
wherein Δz corresponds to the desired voxel size in the z-direction.
The backprojection grid can be arranged in two different ways:
The second method has the advantage that the x,y-positions of the voxels to be reconstructed do not change over λfilt, whereby the backprojection is accelerated; however, an influence (corresponding to the additional interpolation that is required for the inverse Hilbert transformation) can be exerted on the image quality, which must be correspondingly monitored.
In contrast to the standard approach of a filtered backprojection (FBP), the result of a differentiated backprojection (DBP) is not a theoretically exact reconstruction of f(x) but rather instead of this the Hilbert transformation of f(x) along the π-lines described above. The implementation of a two-step Hilbert reconstruction thus contains finding a good way to apply an inverse Hilbert transformation to backprojection results. One method for this is already described by F. Noo et al. in “A two-step Hilbert transform method for 2D image reconstruction”, Phys. Med. Biol. vol. 49, pp. 39093-3923, 2004. In the following the result of the backprojection is designated as (Hf)(x) and the final reconstruction is designated as {circumflex over (f)}(x).
As presented by S. G. Mikhlin in “Integral equations and their applications to certain problems in Mechanics, Mathematical Physics and Technology”, New York: Pergamon, 1957, pp. 126-131, an inverse Hilbert transformation along the direction of a unit vector a can be achieved for a function with limited definition range via the use of
wherein f(x+tα)≡0 for t∉(tmin,tmas). The function C(x) can be calculated in different ways. These are, for example, described by J. Pack et al. in “Cone-beam reconstruction using the backprojection of locally filtered projections”, IEEE Trans. Med. Imag., vol. 24, no. 1, pp. 70-85, January 2005; L. Yu et al., “A rebinning-type backprojection-filtration algorithm for image reconstruction in helical cone-beam CT” in Proc. 2006 IEEE Medical Image Conference (San Diego, Calif.), 2006; and F. Noo et al., “A two-step Hilbert transform method for 2D image reconstruction”, Phys. Med. Biol. vol. 49, pp. 39093-3923, 2004. For our implementation, the simplest way to describe the C(x) function was found with:
since this value can be obtained directly from the measurement data. Equation 12 can now be implemented in Equation 13. The way for this is shown by, for example, J. Pack et al. in “Cone-beam reconstruction using the backprojection of locally filtered projections”, IEEE Trans. Med. Imag., vol. 24, no. 1, pp. 70-85, January 2005 or by F. Noo et al. in “A two-step Hilbert transform method for 2D image reconstruction”, Phys. Med. Biol. vol. 49, pp. 39093-3923, 2004, wherein a rectangular overlap weighting window with a half-pixel shift relative to the initial data is used to avoid aliasing artifacts.
The different variants of the implementation of the differential backprojection in this task can be subdivided into three classes:
Together with the two different methods for generation of the backprojection grids, in total six different versions of the described two-step Hilbert reconstruction method are obtained with differential backprojection and subsequent inverse Hilbert transformation. The different DBP variants used here are subsequently described in detail with regard to their implementation.
applies for a curved detector geometry, with A1(x) and A2(x) as the first and second intersection point of the π-line through (x) with the spiral path and
The z-values of each voxel on the desired s,τ-grid or x,y-grid are first calculated to reconstruct a surface of π-lines. The backprojection for each surface of the π-lines (indexed by Δfilt) is executed over the interval X[λfilt-YFOV/λfilt+n+FOV], wherein YFOV=arcsin(RFOV/RO) with RFOV as the radius of the field of view. A method as it is described by F. Noo et al. in “Exact helical reconstruction using native cone-beam geometries”, Phys. Med. Biol., vol. 48, pp. 3787-3818, November 2003, Chapter 4.3.5, Equations 59 through 64 is used in order to deal with the dependency of the voxels on the region of the projections over which the backprojection is executed. This technique is likewise used in order to determine the limit values in Equation 14.
The implementation of this variant corresponds to that of the Katesevich algorithm that is shown in F. Noo et al., “Exact helical reconstruction using native cone-beam geometries”, Phys. Med. Biol., vol. 48, pp. 3787-3818, November 2003, except for a replacement of the filter step with a derivation according to A given a fixed beam direction.
The third approach to implement the DBP method is based on a rebinning of the measurement data in a pseudo-parallel geometry (“wedge” geometry) corresponding to
wherein ω is initially not changed during the rebinning. The backprojection formula after the rebinning can be described with
and ∂filt=λfilt+π/2 and ∂0=λ0+/2. As is seen here, this approach reduces the filtering to a simpler derivation according to ∂/∂sr, which represents a great advantage for the implementation. Both the backprojection weighting and the voxel dependency of the backprojection range within a given surface of π-lines are eliminated. After deciding on an s,τ- or x,y-backprojection grid and calculation of the z-values of the voxels, the integration can additionally occur directly over the interval [∂filt,∂filt+π].
The control of the CT system ensues with a control and computation unit 10 in which are located computer programs Prg1 through Prgn that can also implement the reconstruction methods according to the invention that are described in the preceding. The output of image data can also ensue via this control and computation unit 10. However, it is noted that the method according to the invention can also be implemented at workstations that merely receive measurement data of a remotely situated CT system.
In summary, the reconstruction method, CT apparatus and computer-readable medium encoded with programming instructions that are disclosed herein implement the image reconstruction along theoretical π-lines, wherein the theoretical π-lines not only lead to interpolated detector data but also can emanate from interpolated source positions. Interpolation thus occurs both at the detector and at the source.
Although modifications and changes may be suggested by those skilled in the art, it is the intention of the inventors to embody within the patent warranted hereon all changes and modifications as reasonably and properly come within the scope of their contribution to the art.
The present application claims the benefit of the filing date of provisional application 60/958,088, filed on Jul. 2, 2007.
Number | Date | Country | |
---|---|---|---|
60958088 | Jul 2007 | US |