System and method for image reconstruction for helical cone beam computed tomography with factorized redundancy weighting

Abstract

A method and system for computed tomography were weighted data is used to reconstruct an image. The apparatus includes an x-ray source producing a cone-beam of x-rays, an x-ray detector disposed to receive x-rays from the x-ray source, a data collection unit, and a processing unit for processing the data using a weighting function using contributions from fan-beam and cone-beam weighting functions. The apparatus can produce images with reduced artifacts.

Description

BACKGROUND OF THE INVENTION

1. Field of the Invention

The present invention relates to computed tomographic (CT) imaging, and in particular to weighting in helical cone-beam CT.

2. Discussion of the Background

As medical CT manufacturers produce scanners with increasing number of detector rows, there arises a need for a practical reconstruction algorithm that can handle the increasing cone angle. Recently an exact helical cone beam algorithm of the shift invariant FBP type (Katsevich algorithm) was proposed, suitable for 1-PI and 3-PI reconstruction. After that practical ways to implement Katsevich algorithm in medical CT scanners were investigated. Generally, exact helical algorithms use only data within the helical PI-intervals, or, equivalently, within the N-PI window [10-11], where N=1, 3, . . . , is the number of helical half-turns used in reconstruction. However from the practical point of view N-PI window weighting has the following properties:

1) Some measured data located outside the N-PI window is not used, which means extra dose to the patient.

2) All data within the N-PI window is used with the same weight; while it is beneficial to the noise reduction, it makes an algorithm more sensitive to patient motion and

imperfections of real data.

3) The N-PI reconstruction restricts the choice of the helical pitch. For example, pitches in the range of 0.75-0.85 are too fast to be used with the 3-PI window, and are suboptimal to use with the 1-PI window, since only a small fraction of data is utilized.

On the other hand, 2D fan beam redundancy weighting has the following advantages:

1) Easily adjusted to the helical pitch

2) Smooth transition from 0 to 1 makes an algorithm more stable to patient motion and imperfections of real data.

SUMMARY OF THE INVENTION

The present invention is directed to a CT method and apparatus where, in one embodiment, the apparatus includes an x-ray source, an x-ray detector disposed to receive x-rays from the x-ray source, a data collection unit, a processing unit for processing the data using a weighting function given as

$w_{\mathrm{CB}}\left(\beta ,\gamma ,v\right)=\frac{u_{\mathrm{FB}}\left(\beta \right)·u_{\mathrm{CB}}\left(v,\gamma \right)}{\sum _{n=-N}^N\left(u_{\mathrm{FB}}\left({\beta }_n\left(\beta ,\gamma \right)\right)·u_{\mathrm{CB}}\left(v_n\right)\right)}$

where β is a projection angle of said x-rays, γ is a fan angle of said x-rays, and ν is a detector coordinate parallel to axis of rotation of said x-ray source.

In another embodiment, the method includes exposing a subject to x-rays from an x-ray source, collecting data, weighting the data using

$w_{\mathrm{CB}}\left(\beta ,\gamma ,v\right)=\frac{u_{\mathrm{FB}}\left(\beta \right)·u_{\mathrm{CB}}\left(v,\gamma \right)}{\sum _{n=-N}^N\left(u_{\mathrm{FB}}\left({\beta }_n\left(\beta ,\gamma \right)\right)·u_{\mathrm{CB}}\left(v_n\right)\right)}$ and reconstructing an image of the subject using the weighted data

BRIEF DESCRIPTION OF THE DRAWINGS

A more complete appreciation of the invention and many of the attendant advantages thereof will be readily obtained as the same becomes better understood by reference to the following detailed description when considered in connection with the accompanying drawings, wherein:

FIG. 1 is a diagram of a system according to the invention;

FIG. 2 is a diagram of scanning geometry;

FIG. 3 illustrates the fan-beam and cone-beam weighting functions;

FIG. 4 is a conventional image of a disk phantom;

FIG. 5 is an image of the disk phantom according to the invention;

FIG. 6 is a conventional image of a chest phantom;

FIG. 7 is an image of the chest phantom according to the invention; and

FIG. 8 is a diagram of conventional filtering and filtering according to the invention.

DETAILED DESCRIPTION

FIG. 1 shows an x-ray computed topographic imaging device according to the present invention. The projection data measurement system constituted by gantry 1 accommodates an x-ray source 3 that generates a cone-beam of x-ray flux approximately cone-shaped, and a two-dimensional array type x-ray detector 5 consisting of a plurality of detector elements arranged in two-dimensional fashion, i.e., a plurality of elements arranged in one dimension stacked in a plurality of rows. X-ray source 3 and two-dimensional array type x-ray detector 5 are installed on a rotating ring 2 in facing opposite sides of a subject, who is laid on a sliding sheet of a bed 6. Two-dimensional array type x-ray detector 5 is mounted on rotating ring 2. Each detector element will correspond with one channel. X-rays from x-ray source 3 are directed on to subject through an x-ray filter 4. X-rays that have passed through the subject are detected as an electrical signal by two-dimensional array type x-ray detector 5.

X-ray controller 8 supplies a trigger signal to high voltage generator 7. High voltage generator 7 applies high voltage to x-ray source 3 with the timing with which the trigger signal is received. This causes x-rays to be emitted from x-ray source 3. Gantry/bed controller 9 synchronously controls the revolution of rotating ring 2 of gantry 1 and the sliding of the sliding sheet of bed 6. System controller 10 constitutes the control center of the entire system and controls x-ray controller 8 and gantry/bed controller 9 and x-rays are emitted continuously or intermittently at fixed angular intervals from x-ray source 3.

The output signal of two-dimensional array type x-ray detector 5 is amplified by a data collection unit 11 for each channel and converted to a digital signal, to produce projection data. The projection data that is output from data collection unit 11 is fed to processing unit 12. Processing unit 12 performs various processing using the projection data. Unit 12 performs data sampling and shifting (described in more detail below), filtering, backprojection and reconstruction, as well as other desired operation on the projection data. Unit 12 determines backprojection data reflecting the x-ray absorption in each voxel. In the circular scanning system using a cone-beam of x-rays as in the first embodiment, the imaging region (effective field of view) is of cylindrical shape of radius R centered on the axis of revolution. Unit 12 defines a plurality of voxels (three-dimensional pixels) in this imaging region, and finds the backprojection data for each voxel. The three-dimensional image data or tomographic image data compiled by using this backprojection data is sent to display device 14, where it is displayed visually as a three-dimensional image or tomographic image.

The system geometry for explaining the apparatus and method according to the invention is shown in FIG. 2. The apparatus and system are directed to a cone-beam system. X-rays are emitted at a plurality of projection angles around the object to collect at least minimally required tomographic data. At each projection angle β, x-rays are emitted with a cone angle γ and a fan angle α and pass through an object with a scanned field of view. The source follows a helical trajectory. The rays are detected by the two-dimensional detector.

A two-dimensional detector can be described where each detector element k, k=1 . . . Nseg×M, where Nseg is the number of detector rows, and M is the number of elements per detector row.

Intensity of the x-ray photon beam (ray) at the detector element k, attenuated by an object or patient, is given by:I_k=I_k⁰ exp(−∫μ(x)dx),  (1)where μ(x) is the attenuation function to be reconstructed, I_k⁰ is the beam intensity before attenuation by μ(x), as produced by the x-ray tube and after penetrating through the x-ray (wedge, bowtie) filter, and ∫μ(x)dx is the line integral of μ(x) along the line l. Mathematically, μ(x) can be reconstructed given a set of line integrals corresponding to a plurality of lines l. Therefore, measured intensity data are to be converted into line integrals first∫μ(x)dx=ln(I_k⁰)−ln(I_k)  (2)

X-ray tomographic reconstruction consists of the following three main steps, data acquisition, data processing and data reconstruction. In data acquisition, the x-ray intensity data are collected at each detector element and each predefined angular view position. This is done within rotating part of the gantry 1. Detector 5 measures incident x-ray flux and converts it into an electric signal. There are two main types of detectors: energy (charge) integrating and photon counting. The electric signal is transferred from a rotating part of the gantry 1 to stationary part though the slipring 2. During this step data may be compressed.

In data processing, data is converted from x-ray intensity measurements to the signal corresponding to line integrals according to equation (2). Also, various corrections steps are applied to reduce effects of undesired physical phenomena, such as scatter, x-ray beam hardening, compensate non-uniform response function of each detector element, and to reduce noise.

Depending on the algorithm, data reconstruction contains all or some of the following processing steps:

• • Cosine (fan angle, cone angle) weighting (can be ×cos, or 1/cos)
  • Data differentiation. It can be performed with respect to fan angle, cone angle, projection angle, source trajectory coordinate, vertical detector coordinate, horizontal detector coordinate, or any combination of those.
  • Data redundancy weighting. Data is multiplied by the weight function W, which may be a function of fan angle, cone angle, projection angle, source trajectory coordinate, vertical detector coordinate, horizontal detector coordinate, or any combination of those.
  • Convolution (filtering). This step utilizes convolution kernel. Some algorithms use ramp-based kernel (H(w)=|w|), some use Hilbert-based kernel (h(t)=1/t, h(t)=1/sin(t), H(w)=i sign(w)). Kernels can be adjusted to the fan beam geometry, scaled, modulated, apodised, modified, or any combination of those. Filtering can be applied along non-horizontal directions, in which case data resampling to the filtering directions and possibly inverse resampling to the original detector grid is required.
  • Backprojection. In this step data is projected back in the image domain. Usually, backprojected data is weighted by a distance factor. Distance factor is inversely proportional to the distance L from the x-ray source position to the reconstructed pixel. Distance factor can be proportional to 1/L or 1/L². Also, some additional data redundancy weighting can be applied here on the pixel-by-pixel basis. Also, usually backprojection step includes obtaining data value corresponding to the ray through the reconstructed pixel by either data interpolation or data extrapolation. This process can be done in a numerous variety of ways.The order in which the above steps are applied depends on a specific reconstruction algorithm.

In the present invention, the processing unit performs data redundancy weighting, termed smooth cone beam weighting, in the following manner. Here, β is the projection angle, γ is the fan angle, ν is the detector coordinate parallel to the axis of rotation, α is the cone angle where ν=Rtan α, g (β, γ, ν) is the cone beam data along the helical source trajectory λ(β)=(R cos β, R sin β, βH/2π), with radius R and pitch H. The physical size of the detector limits detector coordinates to −ν_max<ν<ν_max, and −γ_max<γ<γ_max. In reconstructing an image plane P, normally horizontal planes are reconstructed, so P is given by the equation z=z₀. The values of ν are within the scan range [β_start, β_end], which depends on P.

For each data sample (view, ch, seg), and corresponding ray (β, γ, ν) the following steps are performed:

locate a corresponding image pixel x, which is found by an intersection of the reconstruction plane P and ray (β, γ, ν),

find all measured complementary rays (β_n, γ_n, ν_n), where n=N . . . −N, through x and N is the number of helical half turns. The complementary coordinates are given by:

${\beta }_n\left(\beta ,\gamma \right)=\{\begin{array}{c}\beta +2\gamma +n\pi ,n\mathrm{odd}\\ \beta +2n\pi ,n\mathrm{even},\end{array}{\gamma }_n\left(\gamma \right)=\{\begin{array}{c}-\gamma ,n\mathrm{odd}\\ \gamma ,n\mathrm{even},\end{array}{v}_n\left(\beta ,\gamma ,v\right)=\{\begin{array}{c}\Delta z_{n}R/L^C,n\mathrm{odd}\\ \Delta z_{n}R/L,n\mathrm{even}\end{array}$

where Δz_n=Δβ_nH/2; and ρβ_n=β_S−β_n, where β_S is the view angle corresponding to the image slice position (z₀=β_SH/2;), H is the helical pitch, L=Δz R/ν and L^c=2 R cos γ−L, and

weight the data g(β, γ, ν) depending on the ray position, and normalize by the weighted contributions of all complementary rays.

The cone beam weighting function is given by:

$w_{\mathrm{CB}}\left(\beta ,\gamma ,v\right)=\frac{u_{\mathrm{FB}}\left(\beta \right)·u_{\mathrm{CB}}\left(v,\gamma \right)}{\sum _{n=-N}^N\left(u_{\mathrm{FB}}\left({\beta }_n\left(\beta ,\gamma \right)\right)·u_{\mathrm{CB}}\left(v_n\right)\right)}$ Note that in the summation the index n=0 corresponds to the direct ray: β₀=β, and ν₀=ν.

Examples of functions u_FB(β) (fan beam) and u_CB(ν, γ) (cone beam) are shown in FIG. 3. The functions u_FB(β) and u_CB(ν, γ) are called auxiliary weighting functions. Here Δβ is called the smoothing interval. It can be given as a fixed number of views, as a fixed angular range (for example 10°), or as a percentage of the angular scan length (β_end−β_start). It can be as small as 0% or as large as 50%.

The function u_FB(β) is given by:

$u_{\mathrm{FB}}\left(\beta \right)=\{\begin{array}{cc}0,& \beta <{\beta }_{\mathrm{start}}\\ p\left(\frac{\beta -{\beta }_{\mathrm{start}}}{\mathrm{\Delta \beta }}\right),& {\beta }_{\mathrm{start}}\le \beta <{\beta }_{\mathrm{start}}+\mathrm{\Delta \beta }\\ 1,& {\beta }_{\mathrm{start}}+\mathrm{\Delta \beta }\le \beta \le {\beta }_{\mathrm{end}}-\mathrm{\Delta \beta }\\ p\left(\frac{{\beta }_{\mathrm{end}}-\beta }{\mathrm{\Delta \beta }}\right),& {\beta }_{\mathrm{end}}-\mathrm{\Delta \beta }<\beta \le {\beta }_{\mathrm{end}}\\ 0,& \beta >{\beta }_{\mathrm{end}}\end{array}$

Here the function p( ) can be chosen in various ways. Some examples:

$\begin{array}{cc}\mathrm{Linear}\text{:}& p\left(x\right)=x.\\ \mathrm{Polynominal}\text{:}& p\left(x\right)=3x^2-2x^3.\\ \mathrm{Trigonometric}\text{:}& p\left(x\right)=\frac{1}{2}\left(1-\cos \left(\pi x\right)\right)\mathrm{or}p\left(x\right)={\sin }^2\left(\frac{\pi x}{2}\right)\end{array}$ In general, function p( ) is any function that satisfies: p(0)=0, p(1)=1, and p monotonically increases from 0 to 1.

The function u_CB(ν, γ) can be implemented in various ways. For example it can be given by:

$u_{\mathrm{CB}}\left(v\right)=\{\begin{array}{cc}0,& v<-v_{\max }\\ p\left(\frac{v+v_{\max }}{\Delta v}\right),& -v_{\max }\le v<-v_{\max }+\Delta v\\ 1,& -v_{\max }+\Delta v\le v\le v_{\max }-\Delta v\\ p\left(\frac{v_{\max }-v}{\Delta v}\right)& v_{\max }-\Delta v<v\le v_{\max }\\ 0,& v>v_{\max }\end{array}$ Here Δν is called the smoothing interval. It can be given as a fixed length (for example 3.2 mm, or 3.2 segments), or as a percentage of the detector height 2ν_max. It can be as small as 0% or as large as 50%. The function p( ) can be chosen as discussed above.

Another version, u_CB(ν, γ) can be given by:

$u_{\mathrm{CB}}\left(v,\gamma \right)=\{\begin{array}{cc}0,& v<v_{-}\left(\gamma \right)-\delta v_{-}\left(\gamma \right)\\ p\left(\frac{v-v_{-}\left(\gamma \right)+\delta v_{-}\left(\gamma \right)}{2\delta v_{-}\left(\gamma \right)}\right),& v_{-}\left(\gamma \right)-\delta v_{-}\left(\gamma \right)\le v<-v_{-}\left(\gamma \right)+\delta v_{-}\left(\gamma \right)\\ 1,& v_{-}\left(\gamma \right)+\delta v_{-}\left(\gamma \right)\le v\le v_{+}\left(\gamma \right)-\delta v_{+}\left(\gamma \right)\\ p\left(\frac{v_{+}\left(\gamma \right)+\delta v_{+}\left(\gamma \right)-v}{2\delta v_{+}\left(\gamma \right)}\right),& v_{+}\left(\gamma \right)-\delta v_{+}\left(\gamma \right)<v\le v_{+}\left(\gamma \right)+\delta v_{+}\left(\gamma \right)\\ 0,& v>v_{+}\left(\gamma \right)+\delta v_{+}\left(\gamma \right)\end{array}$ where ν₊(γ) and ν₋(γ) are the boundaries of the n-PI window and are given by:

$v_{+}\left(\gamma ,n\right)=\frac{H\left(2\gamma +n\pi \right)}{4\pi \cos \gamma }$ $v_{-}\left(\gamma ,n\right)=\frac{H\left(2\gamma -n\pi \right)}{4\pi \cos \gamma }$ and n is equal to either 1 or 3, depending on the helical pitch. Here δν₊(γ) and δν₋(γ) are called the smoothing intervals. They can be given as a fixed length (for example 3.2 mm, or 3.2 segments), or as a percentage of the detector height 2ν_max. It can be as small as 0% or as large as 50%. It can also be proportional to:δν₊(γ)=C(ν_max−ν₊(γ))δν₋(γ)=C(ν_max+ν₋(γ)),where C is the proportionality constant, 0<=C<=1.

The function p( ) can be chosen as discussed above, except for it should satisfy an additional condition: p(0.5)=0.5.

Images prepared according to the invention are shown in FIGS. 4-7. FIG. 4 shows a conventional image of a disk phantom obtained using a 320 by 0.5 detector and a helical pitch of 136 mm/rev. The image shown in FIG. 5 was obtained under the same conditions but was processed according to the invention. FIG. 5 is of much better quality.

Similarly, FIG. 6 is a conventional image of a chest phantom while FIG. 7 is an image obtained using the present invention. Again, FIG. 7 is of much better quality where the appearance of many artifacts is reduced

One possible modification of the proposed method is when filtering is applied along non-horizontal filtering directions, as shown in FIG. 8.

The present invention may be implemented in software or in hardware. In particular the operation of the processing unit described above can be carried out as a software program run on a microprocessor or a computer. The software can be stored on a computer-readable medium and loaded into the system.

Numerous other modifications and variations of the present invention are possible in light of the above teachings. It is therefore to be understood that within the scope of the appended claims, the invention may be practiced otherwise than as specifically described herein.

Claims

1. A computed-tomography apparatus, comprising: an x-ray source;an x-ray detector disposed to receive x-rays from said x-ray source;a data collection unit;a processing unit for processing said data using a weighting function given as

2. An apparatus as recited in claim 1, wherein said processing uses said weighting function where:

3. An apparatus as recited in claim 1, wherein said processing unit comprises a filtering sub-unit that applies filtering along detector rows.

4. An apparatus as recited in claim 1, wherein said processing unit comprises a filtering sub-unit that applies filtering in directions other than along detector rows.

5. An apparatus as recited in claim 1, wherein said processing uses said weighting function where: βν+(γ)=C(νmax−ν+(γ))βν−(γ)=C(νmax+ν−(γ))v+(γ, n)=H(2γ+nπ)/4πcosγv−(γ, n)=H(2γ−nπ)/4πcosγC is a proportionality constant, 0<C<1, n equals 1 or 3, and H is a helical pitch.

6. A computed tomography method, comprising: exposing a subject to x-rays from an x-ray source;collecting data;processing said data comprising weighting said data using a weighting function given as

7. A method as recited in claim 6, wherein said processing comprises using said weighting function where:

8. A method as recited in claim 6, wherein said processing comprises filtering along detector rows.

9. A method as recited in claim 6, wherein said processing comprises filtering in directions other than along detector rows.

10. A method as recited in claim 6, wherein said processing comprises using said weighting function where: βν+(γ)=C(νmax−ν+(γ))βν−(γ)=C(νmax+ν−(γ))v+(γ, n)=H(2γ+nπ)/4πcosγv−(γ, n)=H(2γ−nπ)/4πcosγC is a proportionality constant, 0<C<1, n equals 1 or 3, and H is a helical pitch.

11. A non-transitory computer-readable medium containing instructions that may be executed by a computer to perform a method comprising: collecting data associated with exposing a subject to x-rays from an x-ray source;weighting said data using a weighting function given as

12. A medium as recited in claim 11, wherein said method comprises using said weighting function where:

13. A medium as recited in claim 11, wherein said method comprises filtering along detector rows.

14. A medium as recited in claim 11, wherein said method comprises filtering in directions other than along detector rows.

15. A medium as recited in claim 11, wherein said method comprises using said weighting function where: βν+(γ)=C(νmax−ν+(γ))βν−(γ)=C(νmax+ν−(γ))v+(γ, n)=H(2γ+nπ)/4πcosγv−(γ, n)=H(2γ−nπ)/4πcosγC is a proportionality constant, 0<C<1, n equals 1 or 3, and H is a helical pitch.