WO2004051431A2

WO2004051431A2 - 3pi algorithm for spiral ct

Info

Publication number: WO2004051431A2
Application number: PCT/US2003/038375
Authority: WO
Inventors: Alexander Katsevich
Original assignee: Research Foundation Of The University Of Central Florida, Incorporated
Priority date: 2002-12-04
Filing date: 2003-12-04
Publication date: 2004-06-17
Also published as: AU2003300812A8; AU2003300812A1; EP1576446A2; WO2004051431A3; EP1576446A4; JP2006524059A

Abstract

Methods and systems for reconstructing images of moving objects being spirally scanned with two dimensional detectors with a 3PI algorithm. The moving objects can be scanned at a rate of up to approximately three times slower than those of pre-existing systems. In a preferred embodiment, the invention allows for a patient on a table moving through a spiral scanner to be slowed down by a factor of up to three, and still use the same size detector array as those in existing spiral scanning systems.

Description

3 PI ALGORITHM FOR SPIRAL CT

This invention claims the benefit of priority to U.S. Provisional Application, Serial No. 60/430,802 filed, December 4, 2002, and is a Continuation-ln-Part of United States Patent Application Serial No. 10/389,534 filed March 14, 2003 which is a Continuation-ln-Part of Serial No. 10/389,090 filed March Continuation-ln-Part of Serial No. 10/143,160 filed May 10, 2002 now U.S. Patent 6,574,299, which claims the benefit of priority to U.S. Provisional Application 60/312,827 filed August 16, 2001.

FIELD OF INVENTION

This invention relates to computer tomography, and in particular to processes, methods and systems for reconstructing three-dimensional images from the data obtained by spiral scans using 3PI algorithm.

BACKGROUND AND PRIOR ART Over the last thirty years, computer tomography (CT) has gone from image reconstruction based on scanning in a slice-by-slice process to spiral scanning. From the 1970s to 1980s the slice-by-slice scanning was used. In this mode the incremental motions of the patient on the table through the gantry and the gantry rotations were performed one after another. Since the patient was stationary during the gantry rotations, the trajectory of the x-ray source around the patient was circular. Preselected slices through the patient have been reconstructed using the data obtained by such circular scans. From the mid 1980s to present day, spiral type scanning has become the preferred process for data collection in CT. Under spiral scanning a table with the patient continuously moves through the gantry that is continuously rotating about the table. At first, spiral scanning has used one-dimensional detectors, which receive data in one dimension (a single row of detectors). Later, two-dimensional detectors, where multiple rows (two or more rows) of detectors sit next to one another, have been introduced. In CT there have been significant problems for image reconstruction especially for two-dimensional detectors. In what follows the data provided by the two-dimensional detectors will be referred to as cone-beam (CB) data or CB projections.

For three-dimensional (also known as volumetric) image reconstruction from the data provided by a spiral scan with two-dimensional detectors, there are two known groups of algorithms: Exact algorithms and Approximate algorithms, that each have known problems. Under ideal circumstances, exact algorithms can provide a replication of an exact image. Thus, one should expect that exact algorithms would produce images of good quality even under non-ideal (that is, realistic) circumstances. However, exact algorithms can be known to take many hours to provide an image reconstruction, and can take up great amounts of computer power when being used. These algorithms can require keeping considerable amounts of cone beam projections in memory. Additionally, some exact algorithms can require large detector arrays to be operable and can have limits on the size of the patient being scanned.

Approximate algorithms possess a filtered back projection (FBP) structure, so they can produce an image very efficiently and using less computing power than Exact algorithms. However, even under the ideal circumstances they produce an approximate image that may be similar to but still different from the exact image. In particular, Approximate algorithms can create artifacts, which are false features in an image. Under certain circumstances these artifacts could be quite severe.

The first group of algorithms that are theoretically exact and have the shift- invariant FBP structure was disclosed in U. S. Patent Application 10/143,160 filed May 10, 2002, now U.S. Patent 6,574,299, which claims the benefit of priority to U.S. Provisional Application 60/312,827 filed August 16, 2001. A shortcoming of these algorithms is that they operate in the 1PI mode. This means that if the detector array is large in the axial direction, then one has to translate the patient through the gantry very quickly in order to use all of the detector. However, rapid motion of the patient is very impractical for obvious reasons. On the other hand, if the patient moves slowly, only part of the detector is used. This leads to undesirable consequences: part of the x- ray dose is wasted, discretization artifacts are enhanced, noise stability is reduced, etc.

SUMMARY OF THE INVENTION A primary objective of the invention is to provide 3PI algorithms for reconstructing images of objects that have been scanned in a spiral fashion with two- dimensional detectors. For image reconstruction at any given voxel these algorithms require a longer section of the spiral than the 1PI algorithms of U. S. Patent Application 10/143,160 filed May 10, 2002, now U.S. Patent 6,574,299, which is incorporated by reference, which claims the benefit of priority to U.S. Provisional Application 60/312,827 filed August 16, 2001. Consequently, the new algorithms allow to slow the patient down by about a factor of three, but still use the same size detector array.

A first preferred embodiment of the invention uses a five overall step process for reconstructing the image of an object under a spiral scan. In a first step a current CB projection is measured. Next, families of lines are identified on a detector according to a novel algorithm. Next, a computation of derivatives between neighboring projections occurs and is followed by a convolution of the derivatives with a filter along lines from the selected families of line. Next, using the filtered data, the image is updated by performing back projection. Finally, the preceding steps are repeated for each CB projection until an entire object has been scanned. This embodiment works with keeping several (approximately 2-4) CB projections in memory at a time and uses one family of lines.

For the second embodiment, different families of lines can be used in combination with keeping several CB projections in memory. Modifications of these embodiments are possible, that will allow keeping only one CB projection in computer memory at a time. This can be done analogously to what was done in U. S. Patent Application 10/143,160 filed May 10, 2002, which is incorporated by reference, which claims the benefit of priority to U.S. Provisional Application 60/312,827 filed August 16, 2001. Consequently, the new algorithms allow to slow the patient down by about a factor of three, but still use the same size detector array.

Further objects and advantages of this invention will be apparent from the following detailed description of the presently preferred embodiments, which are illustrated schematically in the accompanying drawings.

BRIEF DESCRIPTION OF THE FIGURES

Fig. 1 shows a typical arrangement of a patient on a table that moves within a rotating gantry having an x-ray tube source and a detector array, where cone beam projection data sets are received by the x-ray detector, and an image reconstruction process takes place in a computer with a display for the reconstructed image.

Fig. 2 shows an overview of the basic process steps of the invention. Fig. 3 shows stereographic projection of the spiral onto the detector plane Fig. 4 shows the detector plane with various projections and important lines Fig. 5 shows the boundary curve Fig. 6 shows the continuous illumination case, xl R = (0,0.25,0) Fig. 7 shows the interrupted illumination case, xl R = (-0.5,0,0) Fig. 8 shows the points where various critical events occur Fig. 9 shows the distribution of weights inside the 5IP domain in the case of continuous illumination

Fig. 10 shows the distribution of weights inside the 5IP domains in the case of interrupted illumination

Fig. 1 1 shows the family of filtering lines parallel to L₀

Fig. 12 shows the family of filtering lines tangent to r_±1

Fig. 13 shows how to choose filtering lines depending on the location of x

Fig. 14 shows the family of filtering lines tangent to T^

Fig. 15 shows the filtering lines and the associated constants c_m in different cases:

x F_t (top left panel), x e F₂ (top right panel), x e E₃ and above I", (middle panel),

x G_{ (bottom left panel), x e G₂ (bottom right panel).

Fig. 16 shows possible locations of points s ,s₂,s₃

Fig. 17 shows the filtering lines and the associated constants c_m in different cases:

x F_{ (left panel), e E₂ (right panel).

Fig. 18 is a three substep flow chart for finding families of lines for filtering, which corresponds to step 20 of Fig. 2.

Fig. 19 is a seven substep flow chart for preparing for filtering, which corresponds to step 30 of Fig. 2. Fig. 20 is a seven substep flow chart for filtering, which corresponds to step 40 of Fig.

2.

Fig. 21 is a five substep flow chart for the first part of back-projection, which corresponds to step 50 of Fig. 2. Fig.22 is a three substep flow chart for the second part of back-projection, which corresponds to step 50 of Fig.2.

DESCRIPTION OF THE PREFERRED EMBODIMENTS Before explaining the disclosed embodiments of the present invention in detail it is to be understood that the invention is not limited in its application to the details of the particular arrangements shown since the invention is capable of other embodiments. Also, the terminology used herein is for the purpose of description and not of limitation.

This invention is a continuation-in-part of U. S. Patent Application 10/143,160 filed May 10, 2002, now U.S. Patent 6,574,299, which is incorporated by reference, which claims the benefit of priority to U.S. Provisional Application 60/312,827 filed

August 16, 2001.

The invention will now be described in more detail.

THEORETICAL BACKGROUND First we introduce the necessary notations. Let

C:={y R³ :y_] -Rcos(s),y₂ = Rsm(s),y₃ = s(h/2π),seR}, (l) where h > 0 be a spiral, and U be an open set strictly inside the spiral:

UC{XGR³ :X²+X₂ ² <r²}, 0<r<R, (2)

S² is the unit sphere in R³ , and

D_f(y,β):=[f(y + βt)dt, βeS², (3)

R' :(y-x)-ξ = 0), (4)

that is D_f(y,β) is the CB transform of /. Given (x,ξ)eUx(R³ , 0) , let

S_j = S_j(ξ,ξ-x),j = 1,2,..., denote points of intersection of the plane ϊl(x,ξ) with C .

Also, y(s) := dylds .For β e S² , β^λ denotes the great circle {a e S² :a ■ β = 0} .

Fix any e R³ , where / needs to be computed. In order to compute f(x) we will use a section of the spiral of finite extent, which is to be determined later. For now it will be denoted C(x) . The corresponding parametric interval is denoted I(x) . The main assumption about C(x) is the following.

Property CI. (Completeness condition) Any plane through x intersects C(x) at least at one point.

An important ingredient in the construction of the inversion formula is weight function n(s,x, ),a e β^λ(s,x) . The function n can be understood as follows, x and determine the plane Il(x, ), and the weight n assigned to y(s)≡Yl(x, )nC depends on the location of x . In view of this interpretation we assume n(s,x,a) = n(s,x,-a) . The main assumptions about n are the following.

Property Wl. Normalization condition: _jn(s_],x,a) = \ a.e.on S²; (5)

J

Property W2. There exist finitely many C¹ functions a_k(s,x) e β (s,x),s ≡ I(x) , such that n(s,x,a) is locally constant in a neighborhood

of any (s, a), where el(x) and a e β^λ(s,x),a φ _{k k}(s,x) .

Denote

φ(s,x,θ)

-y(s))ⁿ(^s-,^χ-,^a)-> C-=a(s,θ)eβ (s,x), (6) Under assumptions CI, Wl, and W2 the following general inversion formula is derived in the paper by A. Katsevich "A general scheme for constructing inversion algorithms for cone beam CT," International Journal of Mathematics and Mathematical Sciences, Vol. 21, pp. 1305-1321 (2003) :

(7)

where θ_m [0,π) are the points where φ(s,x,θ) is discontinuous, and c_m(s,x) are

values of the jumps:

c_m(s,x) := lim (φ(s,x,θ_m + ε) - φ(s,x,θ_m - ε)). (8)

Unless n is chosen appropriately, the inversion formula is not necessarily of the FBP type.

3PI LINES AND THEIR PROPERTIES Suppose that the x-ray source is fixed at y(s₀) for some _₀ ε R . Since the

detector array rotates together with the source, the detector plane depends on s₀ and is

denoted DP(s₀) . It is assumed that DP(s₀) is parallel to the axis of the spiral and is

tangent to the cylinder y² + y₂ ² = R² (cf. (1)) at the point opposite to the source. Thus,

the distance between y(s₀) and the detector plane is 2R (see Fig. 3). Introduce

coordinates in the detector plane as follows. Let the d_l -axis be perpendicular to the

axis of the spiral, the d₂ -axis be parallel to it, and the origin coincide with the

projection of y(s₀) . Project stereographically the upper and lower turns of the spiral

onto the detector plane as shown in Fig. 3. This gives the following parametric curves:

p + 2π(j-l)≤s-s₀≤2πj-p,j≥\, or (10)

p + 2πj≤s-s₀≤2π(j + l)-p,j≤-\, (11) where p is determined by the radius of support of the patient: p = 2cos^~i(r/R) (cf.

(2)). These curves are denoted _jJ = ±1,±2,... (see Fig.4). x denotes the projection

of x.

Connect an arbitrary source position y(s₀) to all points y(s) on the spiral, where

2π < s - s₀ < 4π or - 4π < s - s₀ < -2π. (12)

This gives two surfaces, which are denoted S³' (s₀) and S³.^1'1 (s_Q) . The region

bounded by the two surfaces and the cylinder

+x₂ =R² is denoted V^31'1 (s₀) . Let x

be fixed. If s₀ is sufficiently small, then S_y' (s₀) is below x . As s₀ increases, two

cases are possible. In the first one, known as continuous illumination, x enters

V^3P'(s₀) through S™ (s₀) and leaves V^3P'(s₀) through S_" (s₀) . Clearly, the above

procedure yields a unique 3PI interval [b₀(x),t₀(x)] and the corresponding unique 3P1

line Ll" (x) . In the second case, known as interrupted illumination, x enters and

leaves V^3P1 (s₀) several times. More precisely, x intersects each of the surfaces

Sf (s₀) and S_L ^3P' (s₀) exactly three times. Therefore, the above procedure now gives

three 3PI lines L^3PI (x),i - 1,2,3. The corresponding values of the parameter are

denoted b_l(x),t_l(x),i = 1,2,3. Due to the symmetry of the spiral, if x leaves (enters)

V^3PI (s₀) when s₀ - b_t (x) , then x enters (leaves) V^3PI (s₀) when s₀ =t_t(x),i- 1,2,3.

Consider the plane x₃ = 0. The boundary between the two cases is shown in Fig.5.

We see that the curve has no self intersections, so it divides the open disk x² +

< 1 into two regions: N, and X, . In the central one, denoted X_] , there is one 3PI line for

each x . In the exterior one, denoted , , there are three 3PI line for each x .

The following properties can be established. If there is only one 3PI line for a point x , then

xeV³'(_So)^s₀eI^3P,(x):=[b₀,t₀], (13a)

2π<t,-b, <4π, i = 1,2,3. (14a)

If there are three 3PI lines for a point x , then

x V^P'(_So)^s₀eI^3P'(x):=[bA₂A[bA_lA[t₂dΛ (13b)

2π<t,-b,<4π, = 1,2,3. (14b)

In terms of the detector plane equations (13a) and (13b) imply that x is

between T₂ and T_₂ if and only if se I^PI (x) . Therefore, using the analogy with the

1PI case, we can call this region on the detector the 3PI window. In the 1PI case, the parametric interval bounded by the endpoints of the 1PI line of x is called the 1PI

parametric interval of x . Similarly, I ^I (x) is called the 3PI parametric interval of x .

AN AUXILIARY CONSTRUCTION

Suppose first that continuous illumination takes place. Consider the following

two curves on the surface of the unit sphere S². This first curve consists of all unit

vectors a orthogonal to L₀ ^3P! (x) and is denoted by A . The other curve consists of

vectors

_a(j) , _seI™_{(x)t (15)}

and is denoted by T . It is not convenient to represent these curves directly on S² , so they will be shown on the plane using spherical coordinates (θ,θ₂) defined by

S²3 a(s) = (cos θ sin θ₂, sin θ_x sin θ₂, cos θ₂ ), -π < #, ≤π,0 ≤ θ₂ <π. (16) Since both a and - define the same plane, we can restrict <9, to any interval of

length π and "glue" the opposite boundaries using the identification

(£, ) -≡ ((0, + .τ)_mod2;r,7r - 0₂). (17)

A typical situation for xl R = (0,0.25,0) is shown in Fig. 6. It is very convenient to

think about points on S² not only as unit vectors, but also as planes. Each a e S² corresponds to a unique plane through x with normal vector . This correspondence is one-to-one if vectors with opposite orientation are identified.

Together A and T split S² into several domains: D D₂,... . By construction, for a

fixed i , the number of points in C^3P/ (x) nTJ(x,a) is the same for all a e D_t . If Z)

contains k intersection points (IPs), it will be called a ATP domain.

Suppose now that interrupted illumination takes place. In a similar fashion,

define several curves on the surface of S² . The first three curves are obtained by considering unit vectors perpendicular to each of the three 3PI axes. They are denoted A_k,k - 1, 2,3 . The second set of curves is obtained by restricting .s- in (15) to the

intervals [b ,b₂] , [Z>₃J,] , and [t₂ , t₃ ] . These curves are denoted 7|₂ , 7^, , and T₂₃ ,

respectively. A typical situation for xl R = (-0.5,0,0) is shown in Fig. 7.

It is clear that if x, and x₂ are close to each other, they will share similar

diagrams. By this we mean that by a smooth sequence of transformations one diagram can be converted into the other, and in corresponding domains the number of IPs and

their distribution over the subintervals in I ^PI stays the same. An essential change is possible only in a neighborhood of x where a "critical event" occurs: three boundaries

intersect each other at one point on S² . These points can be found numerically and the results are shown in Fig. 8. The smallest distance from any of these points to the center of rotation is r0^ 0.618R . Thus, in situations of interest in medical

applications (r_F0V ≤ 0.5R ), only the following two cases can happen: continuous

illumination as shown in Fig. 6 and interrupted illumination as shown in Fig. 7.

CONSTRUCTION OF THE WEIGHT FUNCTION n

For any s e I ^Pl ' (x) determine the weight function n(s,x,a),a e β (s,x) , according to the following rule:

• In HP domains the only IP gets weight n = 1 ;

• In 3IP domains each IP gets weight n = 1/3 ;

• In 5IP domains two IPs get weight n = 0 , and all the remaining IPs get weight n = 1/3 (see Figs. 9 and 10).

FIRST RECONSTRUCTION ALGORITHM Once the weights are found, we have all the ingredients needed for constructing the algorithm. First of all, for image reconstruction at x we take the 3PI

spiral segment of : C(x) = C ^P! (x) and I(x) = I^3PI (x) (cf. Section 1). As follows

from (6), (7), for each s e I^3PI (x) the filtering directions are determined by finding

the discontinuities of φ(s,x,a) = sgn(a ^■ y(s))n(s,x,a) . The study of these

discontinuities leads us to the following three families of lines.

The first family consists of lines parallel to the spiral tangent. This family is denoted ?₀ . The lines and the associated coefficients c_m are shown in Fig. 1 1 .

Obviously, for any x there is only one line L e £₀ that contains x .

In Fig. 12 we see the second family of filtering lines. It consists of lines tangent to r_±1 and is denoted ?_t . Here for each x within the 1 PI window we take two lines from £ . These lines are determined from the rule explained in Fig. 13. This

fig cure also shows the associated coefficients c r„n .

In Fig. 14 (left panel) we see the family of filtering lines tangent to T . This

family is be denoted £₂. Now, for any x e F u E₂ E₄ E₅ there might be more than

one line Z, e L₂ that contains . The unique line is determined from the following rule

(see Fig. 14, right panel). If x e E, E₄ , the point of tangency is to the right of x . If

e F₂ E₅ , the point of tangency is to the left of . In all cases c_m = 1/3 .

Fig. 15 summarizes the information contained in Figs. 11-14. It shows all possible cases where x might be and all the associated filtering directions and constants c_m . In all cases the direction of filtering as assumed to be from left to right.

It follows from our construction (cf. Fig. 15) that θ_m(s,x) and c_m(s,x) in the general

reconstruction formula (7), (8) depend on x only via β(s,x) . Therefore, we can

replace x by β(s,x) in the arguments of c_m and a¹ and rewrite (7) in the form

1 1 M(s,β) W = -ττ TTT ∑ c_m(s,β) Ψ_m(s,β(s,^χ))ds, (is)

4 4π A ^J *^!"' MW | xJ -- v((sJ) ) | ,

fi" d aγ

V_m(s,β) := [^π —D_f(y(qA sγβ + s\n γA(s,β,θ_m))\ - (19) sin v

x² + x² ≤ (0.618R)². (20)

Step 31. Fix a line L from the said set of lines obtained in Step 20. Note also that given a filtering line L , all x whose projections belong to L and satisfy the rules mentioned above share the same filtering line L (cf. Fig. 15). Thus, (19) becomes a convolution, (18) becomes backprojection, and the algorithm (18), (19) is of the convolution-based FBP type. This can be seen similarly to U. S. Patent Application 10/143,160 filed May 10, 2002, which is incorporated by reference. Other methods and techniques for backprojection can be used. See additionally, for example, US Patent 6,574,299 to Katsevich, which is incorporated by reference.

IMPROVED RECONSTRUCTION ALGORITHM Let ψ(t) be any smooth function defined on R and with the properties

<//(0) = 0 , ψ'(t) > 0J e R . Define the new family of lines £₂ by requesting that any

given line L e £₂ has three points of intersection with T_±1 T^ : s s₂,s₃ , and these

points satisfy: s - s = ψ(s₃ - s₂ ), s + 2π < s < s + 4π, (21) s₃ - s₂ = ψ(s_l -s), s - 4π < s₃ < s - 2π. (22) Recall that _? is the current source position. The lines L e £₂ can be parameterized,

for example, by ->₃,27r <| s₃ \< 4π . Location of _?, and s₂ depends on where s₃ is and

is illustrated in Fig. 16.

Using the properties of ψ it is easy to establish that for each

x ε F_t u E₂ E₄ E₅ there is a unique L e £₂ that contains . Also, if x — >• L" , then

s₂,s₃ -> 2Δ and -?, — > s . Similarly, if x ->^• E_₂ , then s₂,s₃ -> -2Δ and 5, - • s .

More importantly, in the case x e E, E₂ E₄ E₅ it is possible to reduce the

number of filtering directions from two to one. The additional benefit is the improved detector usage. Thus, the top two panels shown in Fig. 15 should be replaced by the diagrams shown in Fig. 17. The case when x appears below L₀ (i.e., x e E₄ E₅ ) can

be obtained from Fig. 17 by symmetry with respect to the origin in the detector plane. The form of the inversion formula (18), (19), remains the same. The first difference is that M(s,β) = 1 (instead of 2) when x e E, E₂ E₄ E₅ . The second

difference is the direction of filtering, which is now determined from (21 ), (22). GENERAL OVERVIEW OF THE RECONSTRUCTION ALGORITHMS

Fig. 2 shows an overview of the basic process steps 10, 20, 30, 40, 50 of the invention. The steps will now be described.

Step 10. Load the current CB (cone beam) projection into computer memory.

Suppose that the mid point of the CB projections currently stored in memory is y(s₀) .

Step 20. Finding families of lines for filtering.

Using either Fig. 15 or, additionally Fig. 17 in case of the improved algorithm, identify the required families of lines. Then, choose a discrete set of lines from each family. Step 30. Preparing for filtering.

Parameterize points on the said lines selected in Step 20 by polar angle. Using interpolation compute the derivative of the CB data (dldq)D_j(y(q), β) | at points

β on the said lines that correspond to discrete values of the polar angle.

Step 40. Filtering. For each line identified in Step 20 convolve the data for that

line computed in Step 30 with filter 1/sin^ .

Step 50. Back-projection. For each reconstruction point x back-project the

filtered data found in Step 40 according to equation (18). Then go to Step 10, unless there are no new CB projections to process or image reconstruction at all the required points x have been completed. Now we describe the algorithm in detail following the five steps 10-50 shown in Fig.

2.

Step 10. Load the current CB (cone beam) projection into computer memory. Suppose that the mid point of the CB projections currently stored in memory is j>(.s₀) . The detector plane corresponding to the x-ray source located at y(s₀) is denoted DP(s₀) .

Step 20. Finding families of lines for filtering.

The set of lines can be selected by the following substeps 21, 22, and 23.

Step 21. From the family of lines £_Q choose an equidistant set of lines that are

parallel to the spiral tangent and that cover the projection of the region of interest onto the detector plane located between T₂ and T_₂ (see Fig. 1 1).

Step 22. From the family of lines £_x choose a discrete set of lines that are

tangent to T, and T_, (see Fig. 12). The extreme left point of tangency on T,

does not have to extend beyond the projection of the region of interest onto the detector plane. Similarly, the extreme right point of tangency on r_, does not

have to extend beyond the projection of the region of interest onto the detector plane. Step 23. From the family of lines £₂ choose a discrete set of lines that are

tangent to T₂ and T_₂ (see Fig. 14, left panel). In both cases the points of

tangency do not have to extend beyond the projection of the region of interest onto the detector plane. In case of the improved algorithm, instead of the lines tangent to T₂ and T_₂ ,

we choose a discrete (say, equidistant) set of values for s₃ on the curves T₂

and r_₂ and then determine the lines L e £₂ by solving equations (21 ), (22).

On both curves the points s₃ do not have to extend beyond the projection of

the region of interest onto the detector plane. Step 30. Preparing for filtering

Step 31. Fix a line L from the said set of lines obtained in Step 20.

Step 32. Parameterize points on the said line by polar angle γ in the plane

through y(s_Q) and L .

Step 33. Choose a discrete set of equidistant values γ that will be used later

for discrete filtering in Step 40. Step 34. For each γ find the unit vector β which points from y(s₀) towards

the point on L that corresponds to γ .

Step 35. Using the CB projection data D_f(y(q),Θ) for a few values of q

close to s₀ find numerically the derivative (dldq)D_j (y(q), Θ) | for all

Step 36. Store the computed values of the derivative in computer memory. Step 37. Repeat Steps 31-36 for all lines L identified in Step 20. This way we will create the processed CB data corresponding to the x-ray source located at y(s₀) .

Step 40. Filtering

Step 41. Fix a line L from one of the families of lines £_m identified in Step

20. Step 42. Compute FFT of the values of the said processed CB data computed in Step 30 along the said line. Step 43. Compute FFT of the filter Msmγ

Step 44. Multiply FFT of the filter 1 /sin γ (the result of Steps 43) and FFT of the values of the said processed CB data (the result of Steps 42). Step 45. Take the inverse FFT of the result of Step 44. Step 46. Store the result of Step 45 in computer memory. Step 47. Repeat Steps 41-46 for all lines in the said families of lines. This will give the filtered CB data Ψ_m(s₀,β ) , where m stands for the line family

number from which L was selected, m - ,\,2 .

By itself the filtering step is well known in the field and can be implemented, for example, as shown and described in U.S. Patent 5,881,123 to Tarn, which is incorporated by reference.

Step 50. Back-projection

Step 51. Fix a reconstruction point x , which represents a point inside the patient where it is required to reconstruct the image.

Step 52. If s₀ belongs to I ^P1 (x) , then the said filtered CB data affects the

image at x and one performs Steps 53-58. If s₀ is not inside the

interval I^3PI (x) , then the said filtered CB data is not used for image reconstruction at . In this case go back to Step 51 and choose another reconstruction point. Step 53. Find the projection x of x onto the detector plane DP(s₀) and the

unit vector β(s₀,x) , which points from y(s_Q) towards x .

Step 54. Using Figs. 11, 13, and the right panel of Fig. 14 identify the lines from the said families of lines and points on the said lines that are close to the said projection . If x is above L₂ or below L"₂ , then in the case of the

improved algorithm use equations (21), (22) and Fig. 16 to find the filtering lines close to the said projection x . This will give a few values of i'_m(s₀,β_l)

for β close to β(s₀, x) .

Step 55. With interpolation estimate the value of Ψ_m(s₀, β(s₀, x)) from the

said values of Ψ„ (J_D, ^ ) for β close to β(s₀, x) .

Step 56. Compute the contribution from the said filtered CB data to the image being reconstructed at the point x by multiplying Ψ_m(s_ϋ, β(s₀, x)) by

-c_m(s,β(s₀,x))/[4π² \ x -y(s₀) |] . The appropriate backprojection coefficient

c_m is selected using Figs. 1 1 , 13, and the right panel of Fig. 14 (see also Fig. 15

for the summary). If x is above If₂ or below L"₂ , then in the case of the

improved algorithm use Fig. 17 to find the appropriate backprojection

coefficient c_m .

Step 57. Add the said contributions to the image being reconstructed at the point x according to a pre-selected scheme (for example, the Trapezoidal scheme) for approximate evaluation of the integral in equation (18). Step 58. Go to Step 51 and choose a different reconstruction point . If all reconstruction points x have been processed, go to Step 59. Step 59. Go to Step 10 and load the next CB projection into computer memory.

The image can be displayed at all reconstruction points x for which the image reconstruction process has been completed (that is, all the subsequent CB projections are not needed for reconstructing the image at those points). Discard from the computer memory all the CB projections that are not needed for image reconstruction at points where the image reconstruction process has not completed. The algorithm concludes when the scan is finished or the image reconstruction process has completed at all the required points.

For example, if the detector does not provide all the data which is required for the 3PI algorithm, one can employ various techniques for estimating the missing data. In this case by the detector (respectively, cone beam projection) we mean the virtual detector (respectively, virtual cone beam projection), which includes both measured and estimated data. If one is able to estimate the missing data exactly, then the 3PI algorithm will produce exact reconstruction. In this sense we still talk about exact reconstruction under real circumstances, when missing data are found approximately.

While the invention has been described, disclosed, illustrated and shown in various terms of certain embodiments or modifications which it has presumed in practice, the scope of the invention is not intended to be, nor should it be deemed to be, limited thereby and such other modifications or embodiments as may be suggested by the teachings herein are particularly reserved especially as they fall within the breadth and scope of the claims here appended.

Claims

I claim:

1. A method of reconstructing images from data provided by at least one detector, comprising the steps of: scanning an object in the spiral fashion with at least one detector that detects at least one cone beam projection, the cone beam projection being wider in the axial direction than projections of four turns of the spiral that are adjacent to a current source position; and reconstructing an exact image of the scanned object in an efficient manner with a convolution based FBP (Filtered Back Projection) algorithm.

2. The method of claim 1, wherein the scanning step includes acquiring two- dimensional cone beam (CB) projection data of the object using the detectors.

3. The method of claim 2, further comprising the step of: using the detectors substantially similar to those required for a 1 PI algorithm.

4. The method of claim 2, wherein the scanning step further includes the step of: detecting the cone beam projection being wider in the axial direction as compared to a cone beam projection used in a 1 PI algorithm.

5. The method of claim 4, wherein the scanning step further includes the step of: detecting the cone beam projection being wider by a factor of at least three times in the axial direction as compared to a cone beam projection used in a 1PI algorithm.

6. The method of claim 1, wherein the object includes: a person.

7. A method of computing exact images derived from spiral computer tomography scan with area detectors, comprising the steps of:

(a) collecting cone beam (CB) projection data from a detector, which is wider than what is required for a 1PI algorithm; the cone beam covering projections of four turns of the spiral that are adjacent to a current source position;

(b) identifying families of lines on a plane LT intersecting the cone beam projection;

(c) preprocessing the CB projection data;

(d) convolution-filtering said preprocessed CB projection data along said lines;

(e) back projecting said filtered data to form a precursor of said image; and

(f) repeating steps a, b, c, d, e, until an exact image of the object is completed.

8. The method of claim 7, wherein the scan includes an x-ray exposure of the object.

9. The method of claim 7, wherein the steps (a) - (f) include: a 3PI algorithm.

10. A method of computing images derived from computer tomography scan with detectors, comprising the steps of:

(a) collecting cone beam (CB) data from a detector during a scan of an object; (b) identifying three families of lines on a plane DP(s) intersecting the cone beam, wherein s is value of the parameter describing the scan path and corresponding to the current source position, and the three families of lines include:

(bi) a first family of lines parallel to y(s), where y(s) is the direction of the scan tangent at the current source position;

(bii) a second family of lines tangent to T, and T_, , where

r, is the projection of the scan turn defined by s < q < s + 2π onto the

plane DP(s) ;

r_, is the projection of the scan turn defined by s- 2π < q < s onto the

plane DP(s) ; q is the parameter along the scan path which describes the point being projected; (biii) a third family of lines tangent to T₂ and T_₂ , where

r₂ is the projection of the scan turn defined by s + 2π < q < s + 4π

onto the plane DP(s) ;

r_₂ is the projection of the scan turn defined by s - 4π < q < s - 2π

onto the plane DP(s) ;

(c) preprocessing and shift invariant filtering said data along said lines of said three families; (d) back projecting said filtered data to form a precursor of said image; and

(e) repeating steps a, b, c, and d until an image of the object is completed.

11. The method of claim 10, wherein the preprocessing includes calculation of the derivative of the CB data with respect to source position.

12. The method of claim 10, wherein the shift invariant filtering includes convolving the said preprocessed data with filter 1/ sin γ .

13. The method of claim 10, wherein back projecting said filtered data from the first family of lines involves multiplying the said filtered data by the coefficient c_m = 2/3 , when the projection of xonto DP(s) is located between L₂ andZT., , where

L^c ₂ is the line parallel to y(s) and tangent to T₂ ;

L ₂ is the line parallel to y(s) and tangent to T_₂ .

14. The method of claim 10, wherein back projecting said filtered data from lines in the first family of lines involves multiplying the said filtered data by the coefficient c_m = 1 /3 , when the projection of xonto DP(s) is located above L₂ or below A₂ .

15. The method of claim 10, wherein back projecting said filtered data from a line in the second family of lines involves multiplying the said filtered data by the coefficient c_m = 2/3 , when the projection of onto DP(s) is located between T, and T_, and

the point where the line is tangent to T, υF, is inside the IPI parametric interval of

x .

16. The method of claim 10, wherein back projecting said filtered data from a line in the second family of lines involves multiplying the said filtered data by the coefficient c_m = -2/3 , when the projection of x onto DP(s) is located between T, and r_, and the point where the line is tangent to T, T_, is outside the IPI

parametric interval of x .

17. The method of claim 10, wherein back projecting said filtered data from a line in the third family of lines involves multiplying the said filtered data by the coefficient c_m = 1 /3.

18. A method of computing images derived from computer tomography scan with detectors, comprising the steps of: (a) collecting cone beam data from a detector during a scan of an object;

(b) identifying three families of lines on a plane DP(s) intersecting the cone

beam, wherein Λ is value of a parameter describing the scan path and corresponding to the current source position, and the three families of lines include: (bi) a first family of lines parallel to y(s) , where

y(s) is the direction of the scan tangent at the current source position;

(bii) a second family of lines tangent to T, and r_, , where

r, is the projection of the scan turn defined by s < q < s + 2π onto the

plane DP(s) ;

r_, is the projection of the scan turn defined by s - 2π < q < s onto the

plane DP(s) ;

-7 is the parameter along the scan path which describes the point being

projected; (biii) a third family of lines on the plane DP(s) that have at least three points

of intersection s_x, s₂,s₃ with r_±1 and T , where T₂ is the projection of the scan turn defined by s + 2π < q < s + 4π

onto the plane DP(s) ;

r_₂ is the projection of the scan turn defined by s - 4π < q < s - 2π

onto the plane DP(s) ;

(c) preprocessing and shift invariant filtering said data along said lines of said three families;

(d) back projecting said filtered data to form a precursor of said image; and

(e) repeating steps a, b, c, and d until an image of the object is completed.

19. The method of claim 18, wherein the points of intersection s_x, s₂, s₃ are

determined according to the following rules: s_x - s = ψ(s₃ - s₂) if s + 2π < s₃ < s + 4π,

s₃ -s₂ = ψ(s_x - s) if s - 4π < s₃ < s - 2π , where ψ(t) is a function with the properties ^(0) = 0 , ψ'(t) > 0,t e R .

20. The method of claim 18, wherein the preprocessing includes calculation of the derivative of the CB data with respect to source position.

21. The method of claim 18, wherein the shift invariant filtering includes convolving the said preprocessed data with filter 1 /sin γ .

22. The method of claim 18, wherein back projecting said filtered data from lines in the first family of lines involves multiplying the said filtered data by the coefficient

c = 2/3 , when the projection of x onto DP(s) is located between L₂ and Z-.₂ .

23. The method of claim 18, wherein back projecting said filtered data from a line in the second family of lines involves multiplying the said filtered data by the coefficient c_m = 2/3 , when the projection of xonto DP(s) is located between F_x and

r_, and the point where the line is tangent to r, uF, is inside the I PT parametric interval of x .

24. The method of claim 18, wherein back projecting said filtered data from a line in the second family of lines involves multiplying the said filtered data by the coefficient c_m = -2/3 , when the projection of x onto DP(s) is located between T_x

and T_, and the point where the line is tangent to T_x T_, is outside the IPI

parametric interval of x .

25. The method of claim 18, wherein back projecting said filtered data from lines in the third family of lines involves multiplying the said filtered data by the coefficient c_m = 2/3 , when the projection of xonto DP(s) is located above L₂ or below L"₂ .

26. A method of identifying a family of lines used for reconstructing images based on a scan of an object in a computer tomography system, comprising the steps of

(i) fixing current source position s , where -y is a parameter describing the scan path;

(ii) picking a plane DP(s) intersecting the cone beam projection; (iii) choosing three points on the scan path, which are described using values of the parameter as s_x,s₂,s₃, such that:

(iiia) I s -s_x \< 2π , (iiib) either 2π <s₂-s <4π and 2π < s₃ - s < 4π or -4π <s₂-s < -2π

and -4π <s₃-s < -2π ;

(iiic) s_x -s = ψ(s₃ -s₂) if 2π <s₃ -s <4πor s₃ -s₂ =ψ(s -s) if

-4π <s₃ -s<-2π , where ψ(t) is a function with ψ(ι) = 0, ψ'(t)> 0,t e R ;

(iv) projecting the three said points onto DP(s) ; and

(v) drawing a line through the said projections.