Image reconstruction from helical partial cone-beam data

A two-dimensional radiation detector (30) receives divergent ray penetrating radiation along conical or pyramidal ray paths which converge at an apex. The radiation may emanate from one or both of an x-ray tube (16) or radionuclei injected into a human subject. As the radiation detector rotates around the subject, it is repeatedly sampled to generate two-dimensional data sets h.sub.t (u,v) at each of a plurality of samplings t. Each two-dimensional data array is weighted (44) and divided into two components. A first component processor (48) convolves a first component of each two-dimensional array with respect to u for each v. A second component processor (50) processes a second component of each two-dimensional array (FIG. 6) including weighting each component with a geometrically dependent weighting function W from a weighting function computer (46). The processed first and second components are combined (52) and backprojected (54).

BACKGROUND OF THE INVENTION 
The present invention relates to the image reconstruction art. It finds 
particular application in conjunction with reconstructing x-ray 
transmission data from CT scanners which move a partial cone-beam 
radiation source along a helical trajectory and will be described with 
particular reference thereto. It is to be appreciated, however, that the 
present application will also find application in reconstructing 
information from other conical or other three-dimensional x-ray sources, 
such as reconstructing conical transmission or emission data in nuclear 
cameras. 
Conventionally, spiral CT scanners include an x-ray source which projects a 
thin slice or beam of radiation. The x-ray source is mounted for 
rotational movement about a subject who is moving along the axis of 
rotation. An arc or ring of radiation detectors receive radiation which 
has traversed the patient. Data from the x-ray detectors represents a 
single spiralling slice through the patient. The data from the detectors 
is reconstructed into a three-dimensional image representation. 
For faster processing, a pair or more of radiation detectors can be 
disposed next to each other. This enables two or more slices of data to be 
collected concurrently. However, like the single slice scanner, only 
intraslice data is used in the reconstruction process. 
One of the difficulties with such prior art scanners is that they place 
major stress on the x-ray generator. When a solid geometric shape of 
x-rays, such as a cone, are generated, the x-rays pass through a 
volumetric region of the subject. In true cone beam reconstruction, 
truncation of the data is not permitted. These x-rays pass along known 
rays, both within traditional planes and at acute angles through several 
planes. The radiation passing through rays at an angle to the central 
plane were previously lost to collimation. By utilizing the radiation 
previously lost in collimation to generate useful diagnostic information, 
the load on the x-ray generator is reduced. 
However, images reconstructed from data collected along divergent beams 
tend to have artifacts. One way of minimizing the divergent ray artifacts 
is to minimize the number of rings, i.e., limit the width of the cone 
beam. Of course, limiting the width of the cone-beam partially defeats the 
original intent. 
Although the additional radiation supplied by the cone-beam is beneficial 
in imaging, it has the detrimental side effect of increasing patient 
dosage. On the other hand, the high dosage enables a volume to be 
constructed with fewer rotations of the cone-beam. 
The present invention provides a new and improved method and apparatus for 
reconstructing volumetric images from cone and other three-dimensional 
x-ray geometries. 
SUMMARY OF THE INVENTION 
In accordance with the present invention, a new and improved method and 
apparatus for reconstructing an image from helical cone-beam data is 
provided. A source generates penetrating radiation about a subject. A 
two-dimensional radiation detector receives the radiation along a 
plurality of divergent rays. The rays are focused at a common origin 
vertex and diverge in two dimensions. The radiation detector is mounted 
and rotated about the subject along at least a helical arc segment of a 
helical path. A sampling means samples the radiation detector at a 
plurality of angular increments along the helical arc segment to generate 
a plurality of two-dimensional views. Each view includes a two-dimensional 
array of data values, where each data value corresponds to one of the 
divergent rays. Data along rows parallel to a tangent to the helix at the 
vertex of the cone are not truncated but data in orthogonal or other 
directions can be truncated. The data collected for each cone is 
referenced with respect to a local coordinate plane. A weighting processor 
weights the data values of each two-dimensional array in accordance with a 
cosine of an angle between a ray corresponding to the data value and a ray 
normal to the radiation detector or local coordinate plane. A first 
component processor convolves a first component of the weighted data 
values one-dimensionally in the direction of the tangent to create a 
convolved first two-dimensional data component array. A second component 
processor processes a second two-dimensional component array of the 
weighted data values. The second processor weights each second component 
value in accordance with characteristics of the helical arc segment and 
positioned along the helical arc segment. The second processor further 
takes a partial derivative of the data values of the second component, and 
multiplies each data value in accordance with a square of a distance from 
the vertex of the divergent rays to the radiation detector corresponding 
to each data value. An adder sums the first and second processed 
components from the first and second component processors. A 
three-dimensional backprojector backprojects the sum from the adder into a 
volumetric image memory. 
In accordance with a more limited aspect of the present invention, the 
second component processor includes a means for computing a line integral 
of the weighted data, a means for computing a cosine weighted first 
derivative of the line integral, a means for computing a product of the 
cosine weighted line integral and a weighting function, a means for 
computing a partial-derivative of the product, a means for 
two-dimensionally backprojecting the partial derivative to generate a 
backprojection, and a means for weighting the backprojection to form the 
processed second two-dimensional component array. 
In accordance with a more limited aspect of the invention, the weighting 
factor is calculated in accordance with characteristics of the helical arc 
segment. 
In accordance with another aspect of the present invention, a method for 
radiographic diagnostic imaging is provided. Penetrating radiation is 
generated. The radiation is received with a two-dimensional radiation 
detector along a plurality of divergent rays. The rays are focused at a 
common origin vertex and diverge in two-dimensions. The radiation detector 
is mounted and rotated about a subject along at least a helical arc 
segment of a helical path. The radiation detector is sampled at a 
plurality of angular increments along the helical arc segment to generate 
a plurality of two-dimensional views. Each view includes a two-dimensional 
array of data values, where each data value corresponds to one of the 
divergent rays. The data values of each two-dimensional array is weighted 
in accordance with a cosine of an angle between a ray corresponding to the 
data value and a ray normal to the local coordinate planes. A first 
component of the weighted data is convolved to create a convolved first 
two-dimensional data component array. A second two-dimensional component 
array of the weighted data is processed, where each second component value 
is weighted in accordance with characteristics of the helical arc segment 
and positioned along the helical arc segment. A partial derivative of the 
data values of the second component is taken. Each data value is 
multiplied in accordance with a square of a distance from the vertex of 
the divergent rays to the radiation detector corresponding to each data 
value. The first and second processed components are then summed. The sum 
is then backprojected into a volumetric image memory. 
One advantage of the present invention is that it provides for faster 
reconstruction of images of volumetric regions. 
Another advantage of the present invention is that it enables volumetric 
images to be reconstructed with less than a full rotation of an x-ray 
source. 
Another advantage of the present invention resides in reduced image 
artifacts. 
Still further advantages of the present invention will become apparent to 
those of ordinary skill in the art upon reading and understanding the 
following detailed description of the preferred embodiments.

DETAILED DESCRIPTION OF THE PREFERRED EMBODIMENTS 
With reference to FIGS. 1a and 1b, a medical diagnostic imaging apparatus 
10 includes an x-ray detector 12 which receives radiation arriving from a 
subject on a subject support 14. In the preferred CT scanner embodiment, 
radiation emanates from an x-ray tube 16 having radiation. The cone-beam, 
as defined herein, can have a circular cross-section, a square or 
rectangular cross-section, a hexagonal cross-section, or the like. 
The radiation detector 12 of the preferred embodiment includes an array of 
detector elements mounted in a configuration which mimics the 
cross-section of the x-ray beam 20. For example, a dozen square detector 
elements 22 may be mounted in a pattern which approximates a circle. When 
the beam is rectangular or square, the detector elements would be 
positioned in a corresponding square or rectangular pattern. As yet 
another alternate embodiment, a plurality of stationary rings of the 
detector elements can be mounted around the subject. Providing a plurality 
of stationary rings of the detector elements adds cost in the detectors, 
but removes the cost associated with rotating the detector assembly. 
Each of the detector assemblies 22 includes a scintillation crystal 24 
disposed toward the radiation source, an array of photodetectors 26 
disposed to view the scintillation crystal, and an array of integrated 
circuits 28 connected with the array of photodetectors 26. Preferably, the 
scintillation crystal, photodetector, and integrated circuit assembly is 
mounted on a common substrate. The scintillation crystal and 
photodetectors are etched or cut to define a larger number of 
light-sensitive elements, e.g., a 16.times.16 array. Optionally, a 
detector collimator 30 is mounted to the scintillation crystals to limit 
received radiation to radiation travelling along rays from an origin of 
the cone-beam collimator. In the CT scanner embodiment, the origin of the 
collimator is selected to match the focal point of the x-ray tube. 
Alternately, the medical diagnostic apparatus may be a nuclear camera. In 
the nuclear camera, the radiation source includes an emission radiation 
source in the form of a radiopharmaceutical that is injected into the 
subject. Because the cone-beam collimator limits received radiation to 
radiation travelling along a conical array of rays, the resultant data is 
again cone-beam data. Further, a transmission line source may be disposed 
opposite the patient from the radiation detector. In a nuclear camera, the 
transmission line source is commonly a radioisotope. 
As is known in the art, an appropriate means or mechanical mechanism is 
provided for rotating the radiation source and the detector around the 
subject and subject support. In a CT camera, the radiation source commonly 
rotates continuously at a relatively high rate of speed. In a nuclear 
camera, the detector and transmission radiation source, if any, commonly 
rotate in incremental steps. 
With continuing reference to FIGS. 1a-1b and further reference to FIGS. 2 
and 3, the radiation detectors are all sampled concurrently to generate a 
current view which is stored in a current view memory or latch 40. The 
memory or latch 40 also stores an indication of the origin of the cone, 
i.e., the location of the x-ray tube and/or the location of the detector 
assembly. In a CT scanner in which the cone is rotating rapidly, the 
detectors are sampled at very short time intervals. In a nuclear camera, 
the output of the detectors is commonly integrated over a dwell duration 
at each angular position. The data from the current view memory 40 is 
conveyed to a partial cone-beam data memory 42 which stores a source cone 
or fan of data h.sub.t (u,v). More specifically, the partial cone beam 
data h.sub.t represents the line integrals or projections onto the plane 
P.sub.0 along the rays .phi. from the origin or x-ray source. A weighting 
processor 44 weights each partial cone data value h.sub.t (u,v) to 
generate a weighted or data view H.sub.t (u,v). In a preferred embodiment, 
the weighting is cosine weighting of the form: 
##EQU1## 
where the weighting term is the cosine of the ray angle and R is the 
radius of the circle that defines the scan helix. 
The object to be reconstructed is represented by a function f from R.sup.3 
to R, where E is the radius of a circle that defines the scan helix, i.e., 
the distance from the z-axis at the origin .OMEGA. to the origin of the 
cone. The trajectory of the x-ray source is described by a curve .phi. 
which is defined mathematically by a function .phi. from an interval 
.LAMBDA. of R to R.sup.3. In the case of helical scanning, such a function 
.phi. is defined by three coordinate equations .phi.(t)=(.LAMBDA.cos 
.omega.t, .LAMBDA.sin .omega.t, .sigma.t) for t a subset of .LAMBDA.. For 
any direction represented by a unit vector .alpha. in R.sup.3, the 
positive half-line originated at .phi.(t) can be represented 
parametrically by {.phi.(t)+r.alpha..vertline.r .epsilon.[0,.infin.]} so 
that h(.alpha., t) is defined by: 
##EQU2## 
is the integral of the function f along the semi-line starting at .phi.(t) 
in the direction .alpha.. A partial cone-beam data at any instant t.sub.0 
may be represented by h(.alpha., t), with .alpha. in the subset of 
R.sup.3. .phi.(t) is the vertex of the partial cone at the particular 
instant t.sub.0. The definition of the function h can be extended for all 
.alpha. in R.sup.3 with .parallel..alpha..parallel..noteq.0. Such an 
extended version of h is denoted hereinafter by g. The inner product of 
two vectors x,y in R.sup.3 is denoted by &lt;x,y&gt;. The Fourier transform of 
the function f is given by: 
##EQU3## 
or in spherical coordinates: 
##EQU4## 
where S denotes the unit sphere R.sup.3. 
The following equations give a relationship between the Fourier transform 
of a function f and that of a function g: 
EQU G(.xi., t)=F(.xi., &lt;.phi.(t), .xi.&gt;) (5), 
where, 
##EQU5## 
is the three-dimensional Fourier transform of g with respect to the first 
variable .alpha.. From the definition itself of the function F, it can be 
concluded that: 
##EQU6## 
If .theta. is a unit vector in R.sup.3, the slice projection theorem which 
produces the relationship between the Fourier transform of f and the 
one-dimensional Fourier transform of its planar integrals R.sub.0 f, it 
can be concluded that: 
##EQU7## 
Consequently, the expression F(.theta., u)-F(-.theta.,-u) is 0 if 
.vertline.u.vertline.&gt;A for A, because the integrals along planes 
perpendicular to .theta. are 0 if the plane does not intersect the support 
of the function f. From this, it can be concluded that: 
##EQU8## 
In condition 1, it is assumed that for any point x in the support of the 
function f, where f(x) is to be reconstructed from the cone-beam data, 
there exists a sub-interval .LAMBDA..sub.x of .LAMBDA. such that any plane 
going through the support of f intersects the subcurve .phi..sub.x of 
.phi. associated with the sub-interval .LAMBDA..sub.x. From this condition 
and by making the change of variables defined by u=&lt;.phi.(t),.theta.&gt;, 
Equation (10) can be expressed as: 
##EQU9## 
where the redundancy weight function M(.theta.,t) satisfies the following 
multiplicity condition which pertains to the number of times that the 
plane going through .phi.(t) and perpendicular to .theta. cuts the curve 
.phi..sub.x : 
##EQU10## 
If the plane perpendicular to a direction .theta. and going through 
.phi.(t) intersects the sub-curve .phi..sub.x at .phi.(s.sub.0), . . . , 
.phi.(s.sub.t.theta.), then the condition of Equation (12) can simply be 
expressed as: 
EQU M(.theta., s.sub.0)+M(.theta., s.sub.1)+. . . +M(.theta., 
s.sub.t.theta.)=1(13). 
The term f(x) can be computed from its Fourier transform using the 
spherical coordinate system as follows: 
##EQU11## 
Substituting Equation (11) into Equation (14), one obtains: 
##EQU12## 
After interchanging the order of integration, one obtains: 
##EQU13## 
Observing that M(-.theta.,t)=M(.theta.,t), after the integration with 
respect to .rho., one has: 
##EQU14## 
Consequently, the computation of f(x) requires: 
EQU .LAMBDA..sub.x ={t.epsilon..LAMBDA..vertline.&lt;x-.phi.(t),.theta.&gt;=0}(18), 
that is, the set of t is such that .phi.(t) is the intersection of the 
plane going through f and perpendicular to the direction .theta.. The 
expression f(x) can also be written: 
##EQU15## 
Observe that g and G are both homogeneous functions with respect to their 
first variable, and that G(.rho..phi.,t)=G(.phi.,t)/.rho..sup.2, Equation 
(19) can be written as: 
##EQU16## 
The inner double integrals of the right hand side of Equation (20) is only 
the three-dimensional inverse Fourier transform of the product of two 
functions, namely: 
EQU G(.xi.,t)-G(-.xi.,t), and .vertline.&lt;.phi.'(t), 
.xi.&gt;.vertline.M(.xi./.parallel..xi..parallel., t). 
Consequently, the integrals are the three-dimensional convolution of their 
inverse Fourier transform. More precisely: 
##EQU17## 
The above convolution backprojection is complex due to the redundancy of 
the weighting function M and is not appropriate for partial cone-beam 
data. 
For the weighting function M constant, and with such constant set equal to 
1 for simplicity, the kernel of the convolution is given by: 
##EQU18## 
The expression .alpha. and .xi. in the coordinate system (T,U,V) with the 
T axis going through .phi.(t), and the U axis parallel to the vector 
.phi.'(t): 
##EQU19## 
i.e., q is the product of the classical ramp filter for the 
two-dimensional reconstruction with the delta function with respect to the 
variable (T,V) weighted with .parallel..phi.'(t).parallel.. Consequently, 
the expression Q(x,t)=(g(.alpha.,t)*q(.alpha.,t)) (x-.phi.(t)) becomes: 
##EQU20## 
Note that the function g is homogeneous with respect to its first 
variable. The convolution expression is a one-dimensional convolution of 
weighted line integrals or projection data along the direction of the 
tangent .phi.'(t) of the trajectory .phi.. The kernel of the convolution 
is the classical ramp kernel which is used in two-dimensional slice 
reconstruction. Furthermore, the convolution of g(-.alpha.,t) is null 
since the object is of compact support. Thus, if M is constant, i.e., 
every plane goes through the point x at which f(x) is to be reconstructed 
cuts a sub-curve .phi..sub.x at the same number of points, then the 
convolution is a one-dimensional convolution in the direction of 
.phi.'(t). Consequently, data can be collected in a partial cone-beam, as 
long as the data are not truncated in the direction of .phi.'(t), for any 
instance t in .LAMBDA..sub.x. 
Although the exact reconstruction formula requires a non-constant weighting 
function M(.theta.,t) which satisfies the redundancy condition stated 
above, based on the simplicity of the convolution expression, the 
weighting function M can be written as the sum of a constant C and a 
non-constant function N, i.e.: 
EQU M(.theta.,t)=C+N(.theta., t) (27). 
With this decomposition, the convolved data is the sum of two terms. The 
first term corresponds to C; and the second term corresponds to the 
weighting function N. Because the first term requires a one-dimensional 
convolution in the direction .phi.'(t), it can be assumed that the partial 
cone-beam data is collected within a rectangle with one side parallel to 
.phi.'(t). 
In condition 2, it is assumed that at any instant t, the data of the 
partial cone beam is not truncated in the direction .phi.'(t) and there 
exists a positive integer K. The integer K has the property that if at 
time t, the data is truncated along a direction .DELTA., then the plane 
defined by the line .DELTA. and the vertex .phi.(t) either (i) intersects 
.phi..sub.x exactly K times or (ii) intersects .phi..sub.x at another 
vertex .phi.(t') for which the data along the line .DELTA.' orthogonal to 
the plane is not truncated. 
Furthermore, it is assumed that along that side, the data is not truncated. 
The redundancy condition on the function M requires N to satisfy the 
following condition: 
##EQU21## 
Observe that if W(.theta.,t).gtoreq.0 and satisfies the condition: 
##EQU22## 
then the function N can be defined by: 
##EQU23## 
W is chosen such that the computation of the convolved data does not use 
the truncated data. Using the truncated data would allow N(.theta.,t)=0. 
Thus, to simplify the reconstruction, the reconstruction formula is chosen 
to be: 
##EQU24## 
where Q.sub.1 is the convolved data corresponding to C which is the 
one-dimensional convolution of the weighted data with the classical ramp 
function in the direction .phi.'(t) and the second term Q.sub.2 is given 
by: 
##EQU25## 
Thus, it follows that: 
##EQU26## 
The right hand side of Equation (33) is only the product of -i/2.pi., with 
the integral of the directional derivative of propagation data along a 
line in the detector plane and orthogonal to the direction .theta.. 
Consequently, the computation of Q.sub.2 does not require the data along 
that line with the truncated data. To insure that such choice is always 
possible, W is selected accordingly. For this reason, appropriate 
conditions are set on the curve .phi. traversed by the x-ray source or 
origin of the cone-beam relative to the axis z. A spiral or a partial 
spiral extending from A to B satisfies the conditions. 
The weighting processor 44 performs the constant weighting component 
weighting of the data while a second factor processor 46 calculates the 
non-constant weighting factor. The weighted data from the processor 44 is 
conveyed to a first or Q.sub.1 processor 48 which calculates the Q.sub.1 
component and to a second or Q.sub.2 processor 50 which calculates the 
Q.sub.2 component with the non-constant weighting function from the 
weighting function processor 46. A processor 52 computes the sum of 
Q.sub.1 and Q.sub.2 and the sum is backprojected by a three-dimensional 
backprojector 54, generally as described by Equations (19)-(21) and as 
described in greater detail below. 
With reference to FIGS. 4, 5, and 6, we start first in the (T,U,V) 
coordinate system. In this coordinate system Q.sub.1 is evaluated as: 
##EQU27## 
To simplify the coordinate transformation, the center of the coordinate 
system (T,U,V) is chosen to be the orthogonal projection C(t) of .phi.(t) 
onto the z-axis. This coordinate system is simply the translation of the 
coordinate system of FIG. 5, with the translation vector D.sub.96, where 
.tau. is the unit vector along the T-axis. In this new coordinate system, 
the coordinates (T,U,V) of the reconstruction point x can be obtained from 
its original coordinates (X,Y,Z) by multiplying the column matrix (X,Y, 
Z-.sigma.t).sup.t with the 3.times.3 rotation matrix A(t) which consists 
of the coordinates of the unit vectors .tau.,.mu.,.nu. along the T,U, and 
V-axes respectively. For .phi.(t)=(Rcos.omega.t, Rsin.omega.t, .sigma.t): 
EQU .tau.=[cos.omega.tsin.omega.t 0] (35a), 
EQU .mu.=[-R.omega.sin.omega.t R.omega.cos.omega.t.sigma.]/(R.sup.2 
.omega..sup.2 +.sigma..sup.2).sup.1/2 and (35b), 
EQU .nu.=[.sigma.sin.omega.t-.sigma.cos.omega.t R.omega.]/(R.sup.2 
.omega..sup.2 +.sigma..sup.2).sup.1/2 (35c). 
It remains then to evaluate the integral expression of Q.sub.1 using a 
sample from the set of line integrals h.sub.t which are estimated from the 
partial cone-beam projection data. Conforming to the choice of the new 
local coordinate system, h.sub.t (u,v) is defined as the integral of the 
function f along the line from .phi.(t) to the point (u,v) in the plane 
defined by C(t), the U and V-axes of the local coordinate system. For P as 
the point x-(0,u,0) as shown in FIG. 3, by the homogeneity property of g, 
one can express g from the corresponding line integrals as follows: 
##EQU28## 
For p as the intersection of the line from .phi.(t) to P with the 
(C(t),U,V) plane, the coordinates of the point p in the coordinate system 
are (0,w(U-u),wV) with w=R/(R-T). Because 
p=.phi.(t)+w(P-.phi.(t))=(R,0,0)+w(T-R,U-u,V), and because w is the 
constant which makes the first coordinate of p null, then: 
EQU p-.phi.(t)=w(P-.phi.(t)) (37a). 
##EQU29## 
Substituting into Equation (36): 
##EQU30## 
Accordingly, the expression of Q.sub.1 becomes: 
##EQU31## 
When v is defined as equal to wV, and after a change of the variable of 
integration is made, Q.sub.1 becomes: 
##EQU32## 
where H.sub.t (u,v) is defined by Equation (1), where 
R=.parallel..phi.'(t).parallel./2.pi..sup.2. In other words Q.sub.1 (x,t) 
is obtained by convolving H.sub.t (u,v) with the classical ramp kernel 
with respect to the first variable u and evaluating it at wU, where U is 
the second coordinate of x in the local coordinate system. H.sub.t (u,v) 
is the projection data h.sub.t (u,v) weighted with the constant divided by 
the distance from .phi.(t) to the point (u,v) on the plane (C(t),U,V). 
Thus, the term Q.sub.1 is essentially a one-dimensional convolution of the 
weighted data with the classical ramp 2D reconstruction filter. By 
contrast, the computation of Q.sub.2 involves the computation of a 
weighting function N based on a knowledge of the number of intersection 
points of a plane with a portion of the helix .phi., the trajectory of the 
vertex of the cones. Accordingly, before calculating Q.sub.2, one first 
establishes the number of intersection points. 
With reference again to FIGS. 1 and 2, and further reference to FIG. 6, the 
equation for the local coordinate planes (cos.omega.t.sub.0, 
sin.omega.t.sub.0, 0) is: 
EQU xcos.omega.t.sub.0 +ysin.omega.t.sub.0 =0 (41). 
A plane (Q) containing .phi.(t.sub.0) and a point x of the support of the 
object intersects .phi. at another point .phi.(t) if and only if the 
intersection I of the line joining the two vertices .phi.(t) and 
.phi.(t.sub.0) with the coordinate plane (P.sub.0) associated with 
.phi.(t.sub.0) belongs to the intersection line D of the plane (Q) with 
the plane (P.sub.0). Note that I belongs to both planes (Q) and (P.sub.0). 
Conversely, when the intersection line D of (Q) with (P.sub.0) contains 
the point I which is the intersection of .phi.(t.sub.0),.phi.(t) with the 
plane (P.sub.0), then .phi.(t) belongs to (Q), i.e., (Q) intersects .phi. 
at .phi.(t). To count the number of intersection points of the plane (Q) 
besides .phi.(t.sub.0), the number of intersections points of the line D 
with the trace of the helix on (P.sub.0) is counted. 
The trace of .phi. on the local coordinate plane (P.sub.0) associated with 
.phi.(t.sub.0) is the set of intersection points I of the line 
.phi.(t.sub.0).phi.(t) with the plane (P.sub.0) as t varies in an interval 
containing t.sub.0. Similarly, the trace of a point x onto the local 
coordinate plane (P.sub.0) is defined by the point of the intersection M 
at the line .phi.(t.sub.0),x with the plane (P.sub.0). 
The trace of .phi. contains two branches of curves which intersect the 
u-axis at infinity, because the line .phi.(t.sub.0).phi.(t) becomes 
tangent to .phi. as t tends to t.sub.0. The coordinates of the 
intersection point I can be written as a point of the line 
.phi.(t.sub.0).phi.(t) as: 
EQU I=(Rcos.omega.t.sub.0 +rR(cos.omega.t-cos.omega.t.sub.0), 
EQU Rsin.omega.t.sub.0 +rR(sin.omega.t-sin.omega.t.sub.0), (42), 
EQU .sigma.t.sub.0 +r.sigma.(t-t.sub.0)) 
where r is a real number. Because the point I belongs to the local 
coordinate plane, the coordinates of I must satisfy the above equation of 
the local coordinate plane. Accordingly, 
r=(1-cos.omega.(t-t.sub.0)).sup.-1. The Equation (42) coordinates of I are 
with respect to the original coordinate system Oxyz. 
With reference to FIG. 3, to obtain the coordinates (u,v) of I with respect 
to the local coordinate system, it suffices to compute the inner product 
of the vector .OMEGA.I with the unit vectors: 
EQU .mu.=[-R.omega.sin.omega.t.sub.0 R.omega.cos.omega.t.sub.0 
.sigma.]/(R.sup.2 .omega..sup.2 +.sigma..sup.2).sup.1/2 and(43a), 
EQU .nu.=[.sigma.sin.omega.t.sub.0 -.sigma.cos.omega.t.sub.0 R.omega.]/(R.sup.2 
.omega..sup.2 +.sigma..sup.2).sup.1/2 (43b). 
More specifically, u=&lt;.OMEGA.I, .mu.&gt; and .nu.=&lt;.OMEGA.I, .nu.&gt;. Because: 
EQU .OMEGA.I=(Rcos .omega.t.sub.0 +rR(cos.omega.t-cos.omega.t.sub.0), 
EQU Rsin.omega.t.sub.0 +rR(sin.omega.t-sin.omega.t.sub.0), (44), 
EQU r.sigma.(t-t.sub.0)) 
it follows that: 
##EQU33## 
as t.fwdarw.t.sub.0, u.fwdarw..+-..infin., and v.fwdarw.0, for t=t.sub.0 
.+-..pi./.omega., Equations (45a) and (45b) become: 
##EQU34## 
The derivatives of u and v with respect to t are: 
##EQU35## 
For .omega. and .sigma.&gt;0, the derivative of u with respect to t is 
negative, and that the derivative of v is positive. Note FIG. 7. This 
analysis of u and v shows that the trace of the helix on a local 
coordinate plane has a graph in the (u,v) space as illustrated in FIG. 8. 
This analysis suggests the following practical criteria to count the 
number of intersection points of a plane going through a point of 
reconstruction with part of a helix. With reference to FIG. 9, [t.sub.1, 
t.sub.2 ] is an interval containing t.sub.0 and I.sub.1, I.sub.2 is the 
trace of .phi.(t.sub.1) and .phi.(t.sub.2) on the local coordinate planes 
(P.sub.0) associated to t.sub.0, respectively. (Q) is a plane containing 
.phi.(t.sub.0) and a point x on the support of the object. M is the trace 
of x on (P.sub.0) and D is the intersection of (Q) with (P.sub.0). 
Looking to the first criteria, (Q) intersects .phi. at .phi.(t) with t in 
[t.sub.1, t.sub.2 ] if and only if D intersects the portion of the branch 
of the trace of the helix from -.infin. to I.sub.1. Practically, this is 
equivalent to the angle between the half lines Mu' and Md' being less than 
the angle between Mu' and MI.sub.1. 
Looking to the second criteria, (Q) intersects .phi. at .phi.(t) with t in 
[t.sub.0, t.sub.2 ] if and only if D intersects the portion of the branch 
of the trace of the helix from I.sub.2 to .infin.. Practically, this is 
equivalent to the angle between the half lines Mu and Md being less than 
the angle between Mu and MI.sub.2. 
Q.sub.2 as defined in Equations (32) and (33) can be further quantified by 
splitting the unit sphere into the union of two disjoint sphere halves, 
namely S/2 and -S/2, the upper and lower half unit spheres. Equation (32) 
then reduces to: 
##EQU36## 
where: 
EQU P(.theta.,t)=(G(.theta., t)-G(-.theta., 
t)).vertline.&lt;.phi.'(t),.theta.&gt;.vertline.N(.theta.,t) (49). 
Substituting -.theta. for .theta. in the second integral of Q.sub.2 and 
because I (-.theta., t )=-I (.theta., t), Equation (48) becomes: 
##EQU37## 
Hence, Equation (32) becomes: 
##EQU38## 
Integrating with respect to the variable .rho. and replacing 
G(.theta.,t)-G(-.theta.,t) by the above-expression, one obtains: 
##EQU39## 
Equation (51) is similar to Equation (11) of the Defrise and Clack 
reference of record, but with a different weighting function N(.theta., 
t). The following notation is used in transposing their results: 
(u,v) is the local coordinates of the intersection of the local coordinate 
plane (P.sub.0) with the line going through .phi.(t) and parallel to the 
direction .alpha.; 
(u.sub.x, v.sub.x) is the local coordinates of the trace of the 
reconstruction point x; 
(s, .psi.) is the polar parameters of the intersection of the local 
coordinate plane (P.sub.0) with the plane perpendicular to the direction 
.theta. and going through x, and .phi.(t); and, 
LH.sub.t (s,.psi.) is the integral of the function H.sub.t along the line 
defined by the polar parameter (s,.psi.), where s is 
the distance from .OMEGA. to the line, and .psi. is the angle between the 
u-axis and the line, see FIG. 10. 
With this notation, Q.sub.2 can be carried out as follows: 
##EQU40## 
may be computed using Grangeat's fundamental relationship by multiplying 
.vertline.&lt;.phi.'(t),.theta.&gt;.vertline.N(.theta.,t) with the weighted 
partial derivatives of the line integrals of H.sub.t : 
##EQU41## 
with H.sub.t being a function defined on the local coordinate plane in a 
fashion similar to the one defined for computing Q.sub.1, except for a 
constant factor: 
##EQU42## 
In other words, Q.sub.2 is obtained by two-dimensionally backprojecting at 
the trace of the reconstruction point, partial derivatives of weighted 
partial derivatives of line integrals of weighted partial cone-beam data. 
With reference to FIG. 11, a coordinate system .OMEGA.qrs is defined with 
.OMEGA.q along .OMEGA..phi.(t), .OMEGA.r parallel to the intersection line 
D of the color coordinate plane with the plane perpendicular to the unit 
vector .theta., and perpendicular to D. In this coordinate system, the 
coordinates of the vector .alpha. are (.rho.,.sigma.,R)/(.rho..sup.2 
+.sigma..sup.2 +R.sup.2).sup.1/2 and those of the intersection point N of 
the line D with the axis .OMEGA.s are (0,s,0). Because the vector .theta. 
is perpendicular to both vectors (0,s,0) and the vector N.phi.(t) with 
coordinates (0,-s,R), it follows that the coordinates of .theta. are 
proportional to the exterior product of the two vectors. Consequently, the 
coordinates of .theta. are (0,R,-s)/(S.sup.2 +R.sup.2).sup.1/2, and hence: 
##EQU43## 
If a and b are the polar and azimuthal angles of the unit vectors .alpha., 
see FIG. 11, then: 
EQU d.alpha.=sin(a)dadb (57). 
Because: 
##EQU44## 
it follows that the Jacobian is: 
##EQU45## 
Hence, Equation (54) becomes: 
##EQU46## 
where h.sub.t (.rho.,.sigma.)=Rg(.alpha.,t). After integrating with 
respect to .sigma., this becomes: 
##EQU47## 
From this, a constructive form of Q.sub.2 can be developed. In the local 
coordinate system in which integration is carried out with respect to 
.theta. over S/2, the vector .theta. can be obtained by multiplying its 
coordinates with respect to (r,s) coordinate system with the rotation 
matrix introduced by the angle (.pi./2-.psi.). Consequently: 
##EQU48## 
With respect to the local coordinate system: 
##EQU49## 
Consequently: 
##EQU50## 
Moreover, if a and b are the polar and the azimuthal angles of the vector 
.theta., then: 
##EQU51## 
and 
EQU b=.psi. (64b). 
Hence, the Jacobian is: 
##EQU52## 
Consequently: 
##EQU53## 
which becomes: 
##EQU54## 
After integrating with respect to s: 
##EQU55## 
With particular reference to FIG. 6, the Q.sub.2 calculation 50 is divided 
into four parts. First, a processor or means 60 computes line integrals of 
the weighted partial cone-beam data. This is followed by a processor or 
means 62 for estimating the partial derivative of the line integrals with 
respect to the variable s, which represents the distance from the origin 
.OMEGA. of the local coordinate system to the lines in accordance with 
Equation (54). The partial derivatives of the line integrals H.sub.t 
(s,.psi.) are weighted by (R.sup.2 +s.sup.2)/(4.pi..sup.2 s.sup.2). This 
weighting is inversely proportional to cos.sup.2 (a), where a denotes the 
divergence angle between the line .phi.(t).OMEGA. and the plane generated 
by the vertex .phi.(t) and the line parameterized by (s,.psi.). For 
simplicity of computation, the parallel beam geometry described by the 
parameters (s,.psi.) on the local coordinate plane is utilized throughout 
the computation of Q.sub.2. 
A second part 46 computes the weight W associated with each line (s,.psi.) 
in order to insure the normalization of the 3D backprojection when 
reconstructing the images. To insure this normalization, the number of 
intersections of a plane with the portion of the helix in which the data 
are three-dimensionally backprojected for reconstruction are computed for 
each plane defined by (s,.psi.). A processor or means 70 identifies the 
affected volume V.sub.t. More specifically, a processor or means 72 
computes the traces A.sub.z and B.sub.z of the end points of the 
backprojection range for each slice, indexed by c, in the affected volume 
V.sub.t. A processor or means 74 computes the weighting function W.sub.t 
for each view t of each slice z in accordance with: 
EQU W.sub.tz (s,.psi.)=N.sub.tz (s,.psi.).vertline.&lt;.phi.'(t), 
(s,.psi.)&gt;.vertline. (68). 
Due to the fact that the truncated data is not used in backprojection, a 
count of the number of intersection vertices with non-truncated data along 
the plane is also counted. 
With reference to FIG. 12 and continuing reference to FIG. 6, without 
counting the current vertex .phi.(t), the number of intersections of a 
plane defined by the line (s,.psi.) is equal to the number of 
intersections of the line (s,.psi.) with the two branches of the trace of 
the portion of the helix in which the data are backprojected to create a 
slice of the reconstructed volume. This number is constant and equal to 1 
when the line (s,.psi.) cuts the line .DELTA. outside the line segment 
defined by the orthogonal projections of .DELTA. of the traces of the two 
extreme vertices of the portion of the helix. Inside this line segment, 
the number is either 0 or 2. That is, the line (s,.psi.) either intersects 
the two branches of the trace or it does not have a common point with the 
trace at all. Second, because the weight is null for lines along which the 
data are truncated, it is not necessary to do any computation of line 
integrals or other related factors of these truncated lines. 
More specifically, it is predetermined into how many slices the data is to 
be reconstructed, the number of pixels per slice, and the like. For each 
view t, a determination is made which slices are affected by view t, K is 
set 76 accordingly to either 1 or 2 for the reconstruction process. Once K 
is set, C from Equation (27) is known as is the weight of Q.sub.1, which 
gives the range of each slice. At 78 a sufficient range [t.sub.1z, 
t.sub.2z ] which satisfies the condition 2, i.e., that the angle between 
the half lines Mu' and Md' is less than the angle between Mu' and MI.sub.1 
to reconstruct slice C in the reconstruction volume V. From the range, the 
weighting for Q.sub.2 is determined. For each view t, a step or means 80 
computes the traces A.sub.z and B.sub.z of .phi.(t.sub.1z) and 
.phi.(t.sub.2z) for each slice z in the affected volume v.sub.t. A step or 
means 82 computes orthogonal projections a.sub.z and b.sub.z of the traces 
A.sub.z and B.sub.z onto the line .DELTA. which makes an angle .psi. with 
the u-axis and goes through the origin .OMEGA.. 
With reference to FIGS. 6 and 12, a step or means 84 computes the number of 
intersections M.sub.tz (s,.psi.) of the portion of the helix between 
.psi.(t.sub.1z) and .psi.(t.sub.2z) with the plane P(s,.psi.) defined by 
.psi.(t) and the line L(s,.psi.), which is orthogonal to the line .DELTA. 
at a distance s from .OMEGA.. Again, the helix can only intersect the 
plane at 1, 2, or 3 points. M.sub.tz (s,.psi.)=1 or 3, if the line L of 
(s,.psi.) cuts [a.sub.z,b.sub.z ]. Otherwise, M.sub.tz (s,.psi.)=2. The 
step or means 86 computes the number M'.sub.tz (s,.psi.) of the vertices 
among the intersections with non-truncated data along the plane 
P(s,.psi.). A step or means 88 computes the number N.sub.tz (s,.psi.) 
which is supplied to the step or means 74 to compute the weighting 
function W.sub.tz (s,.psi.). The number of intersections is calculated in 
accordance with: 
EQU N.sub.tz (s, .psi.)=0 (69a), 
if data along L(s,.psi.) is truncated, and 
##EQU56## 
otherwise. 
With reference again to FIG. 6, a processor 90 determines a partial 
derivative with respect to s and weights it. More specifically, a step or 
means 92 computes I.sub.tz (s,.psi.)=W.sub.tz (s,.psi.)P.sub.t (s,.psi.). 
A step or means 94 determines the Jacobian J.sub.tz (s,.psi.) as follows: 
##EQU57## 
A processor 96 multiplies, at each pixel (u,v) in the local coordinate 
plane, the backprojection of (u,v) from the processor 90 with the square 
of the distance from the vertex .phi.(t) to the pixel. More specifically, 
a step or means 98 adds up 180.degree. of parallel beam backprojections. 
It performs a parallel beam two-dimensional backprojection of J.sub.tz at 
each (u,v) to compute B.sub.tz (u,v). A step or means 100 computes 
Q.sub.2tz (u,v)=(u.sup.2 +v.sup.2 +R.sup.2)B.sub.tz (u,v). These values of 
Q.sub.2 are supplied to the summation means 52 to be combined with the 
corresponding values of Q.sub.1. 
With reference to FIGS. 13 and 14, the three-dimensional backprojector 54 
is a processor which computes the final reconstruction from the processed 
data, such as the above-described convolved data. The backprojector is an 
estimator of: 
##EQU58## 
where Q represents the processed data. Because the data are collected at 
only a finite number of instances, the integral can be approximated by the 
finite sum: 
EQU f(x)=.DELTA.t.SIGMA.Q(x, t.sub.i) (72), 
Further, the continuous volume can be represented by a finite number of 
voxels of size .DELTA.X.DELTA.Y.DELTA.Z, that is x=(1.DELTA.X, m.DELTA.Y, 
n.DELTA.Z)+(X.sub.0, Y.sub.0, Z.sub.0) for some 1=0, 1, . . . , L-1, m=0, 
1, . . . , M-1, n=0, 1, . . . , N-1. (X.sub.0,Y.sub.0, Z.sub.0) is a 
corner of the represented volume. 
When Q is reduced to Q.sub.1, then Q is estimated by appropriately 
weighting the convolved data H.sub.t. Due to the fact that at each 
instance t, data are collected only along a finite number of rays within a 
partial cone-beam, H.sub.t (u,v) must be approximated from a sample 
{H.sub.t (i.DELTA.u+u.sub.0, j.DELTA.v+v.sub.0) .vertline.i=0, 1, . . . , 
I-1, and j=0, 1, . . . , J-1}. It should further be noted that because 
there is data within a partial cone geometry, the backprojection of the 
processed data does not affect the whole represented volume, but only a 
portion of it. Another important point is that the backprojection is 
conveniently carried out in a new coordinate system defined by (C(t),U,V). 
With particular reference to FIG. 13, a three-dimensional image memory 110 
receives the backprojected views. As each view t is received, an affected 
volume determining means or step 112 determines the part of the volume of 
memory 110 which is affected by the processed view t. A voxel generator 
114 sweeps through all of the voxels in the affected part of the volume. A 
coordinate transform processor 116 converts the absolute coordinate system 
(X,Y,Z) of the voxel address generator into the local coordinate system 
(T,U,V). With continuing reference to FIG. 13 and further reference to 
FIG. 14, a weight generating system includes a first weight generator 118 
which generates a first weighting value w which is equal to R/(R-T). A 
two-dimensional input generator 120 which generates local coordinates of a 
voxel given by the voxel generator, the local address being (u,v) and 
generates weighting factors u and w, where u=wU and v=wV. The u and v 
addresses address a two-dimensional memory 122 in which the 
two-dimensional view H.sub.t *RAM is stored. An input interpolator 124 
uses the input data retrieved from the input memory 122 and fractional 
parts of the new coordinates calculated by the two-dimensional input 
generator 120 to estimate the input data corresponding to a ray going 
through the voxel given by the voxel generator. A multiplier 126 
multiplies the weighting factor w by itself to form a weighting factor 
w.sup.2. A multiplier 128 multiplies the interpolated input data H.sub.t 
*RAM(u,v) by the squared weighting function w.sup.2. An adder 130 
retrieves the current values B(l,m,n) retrieved from the volumetric memory 
110 with the product from the multiplier 128. In this manner, the values 
stored in each voxel is updated with each two-dimensional set of processed 
data. Initially, the three-dimensional output memory 110 has each pixel 
value assigned a nominal, non-zero starting value, e.g., 1. In this 
manner, the value at each voxel is iteratively updated until an accurate 
reconstruction is obtained. 
The invention has been described with reference to the preferred 
embodiment. 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.