Title of Invention

METHODS FOR CREATING A PIXEL IMAGE FROM PROJECTIONS

Abstract METHODS FOR CREATING A PIXEL IMAGE FROM PROJECTIONS In the present invention pixel images (116) are created from projections by backprojecting selected projections to produce intermediate images, and performing digital image coordinate transformations (102) and/or resampling on selected intermediate images. The digital image coordinate transformations (102) are chosen to account for view angles of the constituent projections of the intermediate images and for their Fourier characteristics, so that the intermediate images may be accurately represented by sparse samples. The resulting intermediate images are aggregated into subsets (104), and this process is repeated in a recursive manner until sufficient projections and intermediate images have been processed and aggregated to form the pixel image (116). Digital image coordinate transformation can include rotation (102), shearing, stretching, contractions, and the like. Resampling can include up-sampling, down-sampling, and the like.
Full Text This is a continuation-in-part of Provisional Application Serial No. 60/501,350, filed
September 9, 2003, incorporated by reference in its entirety.
This invention relates to tomography, and more particularly, to methods and apparatus
for creating pixel images from projections.
BACKGROUND OF THE INVENTION
Tomographic reconstruction is a well-known technique underlying nearly all of the di-
agnostic imaging modalities including x-ray computed tomography (CT), positron emission
tomography (PET), singly photon emission tomography (SPECT), and certain acquisition
methods for magnetic resonance imaging (MRI). It also'finds application in manufactur-
ing for nondestructive evaluation (NDE), for security scanning, in synthetic aperture radar
(SAR), radio astronomy, geophysics and other areas.
There are several main formats of tomographic data: (i) parallel-beam, in which the
line-integrals are performed along sets of parallel lines; (ii) divergent-beam, in which the
line-Integrals are performed along sets of lines that diverge as a fan or a cone; and(iii) curved,
in which the integrals are performed along sets of curves, such as circles, ellipses, or other
closed or open curves. One problem of tomographic reconstruction is to reconstruct a 2D
or 3D image from a set of its line-integral projections. Another problem of tomographic
reconstruction is to reconstruct a 3D image from a set of its surface-integral projections, that
is, its integrals on a family of surfaces. For example, the 3D Radon transform involves in-
tegrals of the image on a family of 2D planes of various orientations and distances from the
origin. Some of the problems of tomographic reconstruction, and some of the reconstruc-
tion methods, are described in standard references such as F. Natterer, The Mathematics of
Computerized Tomography. Chichester: John Wiley, 1986; F. Natterer and F. Wubbeling,
Mathematical Methods in Image Reconstruction- Philadelphia: Society for Industrial and
Applied Mathematics, 2001; A.C, Kak and M. Soaney, Principles of Computerized Tomo-
graphic Imaging. New York: IEEE Press, 1988; and S.R. Deans, The Radon Transform and
Some of its Applications. New York: Wiley, 1983.

The method of choice for tomographic reconstruction is filtered backprojection (FBP)
or convolution backprojection (CBP), which use an unweighted (in the parallel-beam or
Radon Transform cases)or a weighted (in most other cases) backprojection step. This step is
the computational bottleneck in the technique, with computational requirements scaling as
N3 for an N x N-pixel image in 2D, and at least as NA for an N x. N x N-vaxel image in 3D.
Thus, doubling the image resolution from N to 2N results in roughly an 8-fold (or 16-fold,
in 3D) increase in computation. While computers have become much faster, with the advent
of new technologies capable of collecting ever larger quantities of data in real time (e.g.,
cardiac imaging with multi-row detectors, interventional imaging), and the proliferation of
3D acquisition geometries, there is a growing need for fast reconstruction techniques. Fast
reconstruction can either speed up the image formation process, reduce the cost of a special-
purpose image reconstruction computer, or both.
The dual operation of backprojection is reprojection, which is the process of comput-
ing the projections of an electronically stored image. This process, too, plays a fundamen-
tal role in tomographic reconstruction. A combination of backprojection and reprojection
can also be used to construct fast reconstruction algorithms for the long object problem in
the helical cone-beam geometry, which is key to practical 3D imaging of human subjects.
Furthermore, in various applications it is advantageous or even necessary to use iterative re-
construction algorithms, in which both backprojection and reprojection steps are performed
several times for the reconstruction of a single image. Speeding up the backprojection and
reprojection steps will determine the economic feasibility of such iterative methods.
Several methods have been proposed over the years to speed up reconstruction. For
example, Brandt et al. U.S. Patent No. 5,778,038 describes a method for 2D parallel-beam
tomography using a multilevel decomposition, producing at each stage an image covering
the entire field-of-view, with increasing resolution. Nillson et al. U.S. Patent No. 6,151,377
disclose other hierarchical backprojection methods, While these systems may have ment,
there is still a need for methods and apparatus that produce more accurate images, and offer
more flexibility between accuracy and speed.
Accordingly, one object of this invention is to provide new and improved methods and
apparatus for computed tomography (CT) scanning.
Another object is to provide new and improved methods and apparatus for CT scanning
that produce more accurate images, and offer more flexibility between accuracy and speed.

SUMMARY OF THE INVENTION
These objects are achieved or exceeded by the present invention. Pixel images are cre-
ated from projections by backprojecting selected projections to produce intermediate images,
and performing digital image coordinate transformations and/or resampling on intermediate
images. The digital image coordinate, transformations are chosen to account for view angles
of the constituent projections of the intermediate images and for their Fourier characteris-
tics, so that the intermediate images may be accurately represented by sparse samples. The
resulting intermediate images are aggregated into subsets, and this process is repeated in a
recursive manner until sufficient projections and intermediate images have been processed
and aggregated to form the pixel image.
Digital image coordinate transformation can include rotation, shearing, stretching,
contractions, etc. Reampling can include up-sampling, down-sampling and the like.
Projections can be created from a pixel image by performing digital image coordinate
transformation and/or resampling and/or decimation and re-projecting the final intermediate
image.
BRIEF DESCRIPTION OF THE DRAWINGS
The above mentioned and other features of this invention and the manner of obtaining
them will become more apparent, and the invention itself will be best understood by reference
to the following description of an embodiment of the invention taken in conjunction with the
accompanying drawings, in which:
Fig. 1 is a block diagram of apparatus used for the present invention;
Figs. 2A, 2B and 2C are diagrams of sampling patterns used in some embodiments of
the present invention;
Figs. 3A, 3B and 3C are additional sampling patterns used in some embodiments of
the present invention;
Fig. 4 is a diagram illustrating a known method of backprojection;
Fig. 5A is a diagram illustrating an algorithm for one embodiment of the present
invention:

Fig. 5B is a diagram illustrating the manner in which intermediate images axe gener-
ated in the embodiment of Fig. 5A;
Fig. 6 is a diagram illustrating Fourier characteristics used to produce the intermediate
images of Fig. 5A;
Figs. 7A, 7B and 7C are diagrams showing Fourier supports of intermediate images
for the backprojection algorithm illustrated in Fig. 5A, when the coordinate transformation
is a digital image rotation;
Fig. 8 is a diagram illustrating an algorithm used in another embodiment of the present
invention;
Fig. 9 is a diagram showing the evolution of the spectra] support in the algorithm of
Fig. 8, the blocks (1...9) corresponding to the corresponding points in Fig. 8;
Fig. 10A is a diagram describing an algorithm for embodiment of the present inven-
tion;
Fig. 10B illustrates the image coordinate transformation used in the embodiment of
Fig. 10A:
Fig. 1 ] A is a diagram illustrating shear scale backprojection, and Fig. 11B is a dia-
gram illustrating hierarchical shear scale backprojection;
Figs. 12A and 12B are diagrams showing the effect of image shearing on the spectral
support of intermediate images;
Fig. 13 is a diagram illustrating an algorithm for another embodiment of uie present
invention;
Fig. 14 is an algorithm for finding optimal shear factors;
Fig. 15 illustrates an algorithm for another embodiment of the present invention;
Fig. 16 illustrates an algorithm for still another embodiment of the present invention;
Fig. 17 is a diagram which illustrates common fan beam geometry with a circular
scanning trajectory;
Fig. 18 illustrates an algorithm for another embodiment of the present invention;
Fig. 19 illustrates weighting of intermediate images in the algorithm described in Fig.

38;
Figs. 20A, 20B and 20C illustrate sampling points for the second hierarchical level of
Fig. 18;
Figs. 21A and 21B arc diagrams of sampling patterns used in the algorichm of Fig. 18:
Fig. 22 is a diagram showing original intersection points obtained using the method
illustrated in Figs. 20A-20C;
Figs. 23A, 23B, 23C and 23D illustrate sampling points for rotation for up-sampling
used in Fig. 18;
Fig. 24 illustrates local spectral support at a point of an intermediate image generated
by the algorithm of Fig. 18;
Fig. 25 is a diagram of nonuniform sampling patterns used in the algorithm of Fig. 18;
Figs. 26A and 26B illustrate sampling patterns for resampling and coordinare trans-
formation;
Figs. 27A and 27B illustrate alternative sampling schemes which can be used in the
present invention;
Fig. 28 includes two diagrams of divergent beams used in the present invention:
Fig. 29A is a diagram showing a conebeam, and Fig. 29B is a diagram illustrating
resampling;
Fig. 30 is a diagram of an algorithm used for resampling projections;
Fig. 31 is another algorithm used for resampling in the embodiment of the present
invention;
Fig. 32 is a diagram of diagram of an algorithm used for fast hierarchical reprojection;
Fig. 33 is a diagram of another algorithm for fast hierarchical reprojection;
Fig. 34 is a graph showing the results of experiments using the present invention;
Fig. 35 includes sample images generated with the present invention;
Fig. 36A is a display of a reconstructed image ' sing the conventional algorithm, and
Fig. 36B shows a result obtained with the fast algorithms of the present invention; and

Figs. 37A and 37B are diagrams of point spread functions of algorithms comparable
to conventional algorithms, and Fig. 37C displays a point spread function of a fast algorithm
of the present invention.
DETAILED DESCRIFTION
Symbols and Fonts
The following system of mathematical symbols and fonts will be used to improve
clarity.
Functions in the space domain are denoted by small letters (e.g. f(x)), while their
Fourier transforms are denoted by capital letters (F(W)).
The indices of two-variable functions are denoted variously, depending on conve-
nience. The following three notations of function f are equivalent: f(x1,x2), f (f), and
f(£]).
Continuous-domain and discrete-domain functions respectively are distinguished by
the style of parentheses used with their indices: f(x1 x2) is a function of two continuous
variables (i.e., f e R2), and f(m1, m2] is the sampled version of f(x) and is therefore a 2-D
array.
A linear operator and its corresponding matrix are distinguished by the font style.
Suppose (A E R2x2) is a matrix., then A is its associated linear operator. For example, if A
is a coordinate transformation, g(f) — (Af) (r) = f{Ax).
The same operator is sometimes denoted differently inside and outside block diagrams.
While outside it may be denoted as A{a), within the block diagram it is denoted as Aa.
Overview of Hardware

The presenc invention has application in a variety of imaging apparatus, including CT
scanners. Typical imaging apparatus 10 (Fig, 11) includes a scanner 12 which acquires data
from an object such as a head, and sends raw daca corresponding to line-integral projections,
e.g., with a divergent beam geometry, 14 to a projection pre-processor 16. The projection
pre-processor 16 applies various conversions, normalizations, and corrections to the data, as
well as weighting and filtering, which may be shift varying. The output of the projection
pre-proccssor 16 is a collection of pre-processed projections, hereinafter simply referred to
as projections, also called sinogram). 18, which is fed to a sinogram update processor 20. The
sinogram update processor 20 possibly modifies the input sinogram1 18, using information
from sinogram2 34, for example correcting for various artifacts including beam-hardening,
or as part of a multi-step or iterative reconstruction procedure.
The output of the sinogram update processor 20 is a sinogram3 22, which , input
to a fast, backprojection processor 24. The fast backprojection processor' 24 is generally
a computer or special purpose hardware, or any combination thereof, of any suitable type
programmed and/or wired to perform the algorithms described herein.
The output of the fast backprojection processor 24 is an electronic image1 26, which
is input to an image conditioning processor 28. The image conditioning processor 28 per-
forms necessary postprocessing of the electronic image, possibly including identification
and extraction of artifact images, or images for further processing in a multi-step or iterative
reconstruction process.
If desired, the image conditioning processor 28 can produce au electronic image, 30
that is fed to a fast reprojection processor 32. The fast reprojection processor 32 is generally
a computet or special purpose hardware, or any combination thereof, of any suitable type
programmed and/or wired to perform the algorithms described herein. If desired, this pro-
cessor can share the use of the same computer and hardware employed by the backprojection
processor 24.
The output of the fast reprojection processor 32 is a sinograms2 34, which is fed back
into the sinogram update processor 20. The backprojection/reprojection process can continue
until suitable results are. obtained. While reprojection is not always needed, it is helpful in
many situations.
When the electronic image1 26 is suitable, the image conditioning processor 28 pro-
duces an electronic image3 36, which is fed to a storage/analysis/display device 38. It is

contemplated that the electronic image 3 36 can be stored in a computer memory, and/or
analyzed electronically for anomalies or dangerous materials, for example, and/or displayed,
and/or printed in some viewable form.
Overview of Backprojection and Reprojection Methods of the Present Invention
The backprojection methods of the present invention use various techniques to create.
an image made of pixels (picture elements) and/or voxels (3D picture elements), hereinafter
referred to collectively as pixels, which will now be introduced in a general way.
This explanation uses terminology and processes commonly used in multi-dimensional
signal processing, for example as described in D. Dudgeon and R, Mersereau, Multidimen-
sional Digital Signal Processing. Englewood Cliffs: Prentice-Hall, 1983. Some terms in this
description of the present invention are used in the following contexts. The term Sampling
pattern refers to a set of points in space with positions defined relative to a system of coor-
dinates. Examples of sampling patterns are seen in Figs. 2A-2C and 3A-3C, A Cartesian
sampling pattern refers to a set of points formed by the intersection of two mutually per-
pendicular sets of parallel lines. The term continuous image refers to a function defined on
, a coordinate system, for example, f(x: y), and f(x, y, z) are respectively 2D and 3D func-
tions. A digital, image is an array of values of a continuous image on a sampling pattern.
More broadly, a continuous image can be represented by an array of numbers that serve
as the coefficients in a series expansion with respect to some basis set, such as splines, of
which the Cartesian product of zero-th order splines yields the familiar square pixel form for
displaying digital images as continuous images. Hereinafter, this array of numbers will be
also referred to as a digital image. All images stored in a digital computing device must be
digital. For brevity, both digital and continuous images will be often referred to hereinafter
simply as images, with the meaning inferred from the context. With this terminology, a pixel
image is a digital image corresponding to a sampling pattern that is a lattice, i.e., a uniformly
spaced, periodic pattern, usually but not necessarily Cartesian.
One sampling pattern will be said to be sparser than another, if it yields a smaller total
number of samples of a continuous image. Typically a sparser sampling pattern will have
a lower sampling density. Oversatnpling refers to using more samples than necessary to

represent a continuous image at a desired accuracy. The corresponding digital image will be
said to be oversampled. Given a digital image corresponding to a continuous image for one
sampling pattern, the process of producing a new digital image corresponding, with a desired
accuracy, to the same continuous image on a different sampling pattern, will be called digital
image resampling. Up-sampling and down-sampling are special cases of resampling on a
denser or sparser sampling pattern, respectively. Further, up-sampling or down-sampling by
a factor of 1 involves no change in the digital image, is considered a form of up-sampling
or down-sampling. Digital image addition refers to to a point-by-point addition of digital
images defined with respect to the same sampling pattern, or the same basis, in the case of
an expansion with respect to a basis. Lower dimensional digital filtering refers to digital
filtering of a multidimensional array along only a subset of its dimensions, for example,
separate ID filtering of each column of a 2D rectangular digital image.
Coordinate transformation of a continuous image refers to operations such as rotation,
shearing, stretching, or contraction. To define a digital coordinate transformation, consider
two continuous images related by a coordinate transformation, and the corresponding digi-
tal images representing the continuous images with respect to a common sampling pattern.
The process of producing one digital image from the o.her is called digital image coordinate
transformation. This can be accomplished by digital filtering, i.e., by discrete index opera
tions on the digital image. Specific examples include (but are not limited to) digital image
rotation, digital image shearing, digital image stretching, digital image contraction, and the combinations of such operations. Methods for performing digital image coordinate transfor-
mation are known, for example, as described in M. Unser, P. Thevenaz, and L. Yaroslavsky,
Convolution-based interpolation for fast, high-quality rotation of images. IEEE Transactions
Image Processing, Vol. 4, pp. 1371-1381, 1995.
Some digital image coordinate transformations are illustrated in Figs. 2A-2C and 3A-
3C. Fig. 2A shows the outline of a continuous image (a rectangle) and the sampling pattern
for a digital image representing it. Values of the continuous image on the heavy dots are
included in the digital image. Figs. 2-B and 2-C show the rotated and the sheared continuous
image, respectively, on the same sampling pattern, With the heavy dots showing the values
included in the digitally rotated/sheared version of the digital image in Fig. 2-A.
Fig. 3A also shows a continuous image and a sampling pattern defining a digital
image. Fig. 3B shows the stretched continuous image by some constant factors in the x and
y dimensions. The digitally stretched image is defined by values of the stretched continuous

image on the heavy dots. Fig. 3C shows the same continuous image as in Fig. 3A. but with
a sampling partern denser by certain stretch factors. The corresponding digital image in Fig.
3C will be the same as in Fig. 3B. More generally, digital image stretching or contraction
can be equivalent, for sampling patterns with some regularity, to digital image up-sampling
or down-sampling.
Note that the application of certain coordinate transformations, such as rotation by 0
degrees, shearing by a shear parameter of zero, or stretching or contraction by a factor of 1,
leave the digital image unchanged, and therefore may be cither included or omitted in the
description of a process, without changing the result.
The Fourier transform of a continuous image allows one to determine, via sampling
theory, the properties of sampling patterns such that the corresponding digital image repre-
sents the continuous image to desired accuracy, as explained, for example, in D. Dudgeon
and R. Mersereau, Multidimensional Digital. Signal Processing. Englewood Cliffs: Prentice-
Hall, 1983. Likewise, the discrete-time Fourier transform (DTFT) of a digital image allows
one to determine what digital image coordinate transformations will produce a digital image
that represents the transformed continuous image to a desired level of accuracy. The relevant-
properties of the Fourier transform of continuous images, and the DTFT of digital images,
will be collectively referred to as Fourier characteristics.
Weighted backprojection operations require weighting of each projection by a weight that depends on the position of the pixel generated by the backprojection. Different weights
can be used for different projections, depending, for example, on the position of the source
at which the projection was acquired, as explained in A.C. Kak and M. Slaney, Principles
of Computerized Tomographic Imaging. New York: IEEE Press, 19S8. As a special case,
the weighting factor can be equal to 1, which corresponds to no weighting, but is a weight-
ing factor. Unless specifically indicated, weighted and un-weighted backprojection will be
collectively referred to as backprojection.
With this background information, several embodiments of the invention will be de-
scribed.
The backprojection processor 2S (Fig. 1) is programmed to perform the algorithms
used to practice the present invention. The algorithms will be discussed in detail, but will
first be described in more general terms. Steps in the backprojection process are indicated in
block diagrams, with reference numerals.

Fig. 4 illustrates rotation-based backprojection, as the sum of rotated, intermediate im-
ages formed by backprojecting individual projections q1~qp at zero angle in step 50. The
backprojections produces images which are subjected to coordinate transformation at 52.
and aggregated at 54 to produce an image /. By itself, this structure is equivalent to con-
ventional backprojection, and offers no reduction in operation count. However, it serves as
a stepping stone to the introduction of some of the fast hierarchical backprojection methods
of the present invention-
Fig. 5A illustrates a hierarchical backprojection method for creating a pixel image
/ from a plurality of projections q1.-.qp. Each projection qm is backprojected at 100 to
produce a plurality of intermediate images I1,1---I1,P.. This is the zeroth or preparatory level
of a hierarchical backprojection. Digital image coordinate transformations are performed
on selected intermediate images at 102. Subsets of the transformed intermediate images
(pairs, here), are aggregated at 104 to produce aggregate intermediate images I2l, ...I2,p/2-
This is the. first level of a hierarchical backprojection. The aggregate intermediate images
of the first level serve as new intermediate images input to the next level of hierarchical
backprojection. The process of applying digital image coordinate transformations to selected
intermediate images at 106, and aggregating selected intermediate images at 108 to produce
new intermediate images continues until all intermediate images have been processed and
aggregated into the final image f at 116.
If desired, the operations within and across some of the levels can be combined. For
example, the zeroth and first level can be combined for some of the projections, and the
corresponding intermediate images from the set I21 > • ■ -I2p/2 produced instead by a backpro-
jection 112 at selected view angles of two or more (exactly two, for the embodiment in Fig.
5A) selected projections qp, as shown in Fig. 5B. Alternatively, some of the initial interme-
diate images can be produced by an equivalent process involving no explicit backprojection,
such as a direct Fourier reconstruction method, as described in F. Natterer and F. Wubbeling,
Mathematical Methods and image Reconstruction. Philadelphia: Society for Industrial and
Applied Mathematics, 2001.
The parameters of the various digital image coordinate transformations are chosen to
account for the view angles of the constituent projections of the intermediate images, and for
the Fourier characteristics of the intermediate images, so that the aggregates of the interme-
diate images may be accurately represented by sparse samples, as will be explained. These
Fourier characteristics focus on the essential spectral support of the intermediate images, i.e.

the region in the Fourier (frequency) domain, where the Fourier transform of the intermediate
image is significantly different from zero. Sampling theory teaches that the spectral support
of a continuous image determines the nature of the sampling patterns that produce digital im-
ages that represent the underlying continuous image, and from which the continuous image
can be reliably reconstructed. In particular, Fig. 6 shows the typical wedge-shaped spectral
support of intermediate image I1m in the hierarchical algorithm.
Figs. 7(A) and 7(B) show Fourier supports of intermediate images in the binary hierar-
chical backprojection algorithm illustrated in Fig. 5(A), when the coordinate transformation
is chosen to be a digital image rotation. Fig. 7(A) shows the Fourier-domain support of
the virtual intermediate images I. A virtual image I1mis composed of the backprojection of
the projections that are included in the corresponding image I1m. Fig. 7(B) shows that by
choosing the parameters of the coordinate transformation to account for the view angles of
the constituent projections of the intermediate images, and for the Fourier characteristics of
the intermediate images, the vertical bandwidth (the height of the broken-line rectangle) of
the intermediate images can be minimized, allowing sparse sampling and reducing compu-
tational requirements. Fig. 7(C) shows the space-domain sampling scheme of I1m,It[_12rn.
and Ig_1i2,„_7. The sampling points are at the intersections of the dotted lines in the space
domain.
Fig. S describes another embodiment of the present invention, which performs a
ternary aggregation of intermediate images, with P = 2 x 3L projections. For concrete-
ness of illustration, the case of P = 18 projections is shown. As in Fig, 5(A), projec-
tions qe1. qe, ".q08 are backprojected at 100 to produce a plurality of intermediate images
I1, ...if ]S. This is the. zeroth level of a hierarchical backprojection.
The projections are grouped in 3's, and selected projections (two of the three in Fig.
8) in each group are subjected to digital image coordinate transformation at 102. Subsets
(triplets, in Fig. 8) of the transformed intermediate images If, are aggregated at 104 to
produce aggregate intermediate images I2,1 ...,I26. This is the first level of hierarchical
backprojection (I = 1).
In the second level of hierarchical backprojection (I = 2) selected aggregate interme-
diate images Id2 m undergo a coordinate transformation composed of stretching or upsampling
along the y coordinate at 106, and another coordinate transformation /C2m on selected ag-
gregate intermediate images at 108. Here, also, the transformed intermediate images are.

aggregated in groups of three at 110 to produce intermediate images if^, I^2. Selected inter-
mediate images from level two are subjected to digital image coordinate transformations al
112 and aggregated at 114 to produce the next level intermediate images I%x. This process is
repeated to the extent necessary to produce the image / at 116, depending on the number of
projections. The digital image coordinate transformation denoted by K, at the last level 112 is
a digital image rotation, and those at the preceding levels 102 and 108 can also be chosen to
be digital image rotations. Here too, the parameters of the various digital image coordinate
transformation are chosen to account for the view angles of the constituent projections of the
intermediate images, and for the Fourier characteristics of the intermediate images, so that
the aggregates of the intermediate images may be accurately represented by sparse samples,
as will be explained.
Fig. 9 shows the evolution of the spectral support in me ternary backprojection al-
gorithm of Fig. 8. The numbers (1),...,(9) correspond to the corresponding points in the
block diagram of the algorithm in Fig. 8. The spectra shown are the discrete-time Fourier
transforms (DTFTs) of the digital images Ilm.
Fig. 10(A) describes an algorithm for another embodiment of the present invention,
using ternary two-shear-based hierarchical backprojection. The embodiment of Fig. 10(A)
is similar to the embodiment of Fig. 8, and the reference numerals from Fig. 8 are used
where appropriate. As in Fig. 8, P = 2 x 3L projections, and the case of L — 2 is shown.
The embodiment of Fig. 10(a) differs from that described in Fig. S in two respects. First, an
additional upsarnpling slep along the x coordinate is included at 101 in the first level coordi-
nate transformations. Second, at least some of the digital image coordinate transformations,
denoted by AT at 102 and 108, are. composed of a sequence of two digital image shear op-
erations, the first (120) along the y coordinate (122), the second along the x coordinate, as
shown in Fig. 10(b).
Fig. HA illustrates shear-scale backprojection, and Fig. 1 IB illustrates an equivalent
hierarchical shear-scale backprojection. In Fig. 11A, the plurality of projections q1:..., q4,
are backprojected at 140 and subjected to a shear-scale coordinate transformation at 142.
The resulting intermediate images are aggregated at 144.
In Fig. 11B, the projections q1:..., q4 are backprojected at 146, subjected to a shear-
scale coordinate transformation at 148, and aggregated in subsets at 150. The intermediate
images produced by the aggregation are subjected tc a shear-scale transformation again at

152, and the resulting intermediate images are aggregated at 154. This process continues in
a recursive manner until all of the projections and intermediate images have been processed
and aggregated to form an image /. Here too, the parameters of the digital shear transforma-
tions and sampling of intermediate images are selected to account for the view angles and the
Fourier characteristics, so as to reduce the sampling requirements and required computation.
Fig. 12 shows the effect of image shearing on the spectral support of intermediate
images. Fig. 12A shows the spectral support of a certain subset of the projections as they
should appear in the final image. Fig. 12(B) shows the spectral support of the intermediate
image composed of the same projections, with coordinate transformation parameters chosen
to minimize the highest radian frequency in the vertical direction.
Fig. 13 illustrates an algorithm for another embodiment of the present invention, ;i
ternary hierarchical shear-scale backprojection. Projections qe1, ...,q0l are backprojected at
160 and subjected to a shear-scale digital image coordinate transformation at 162. Subsets
of the resulting images are aggregated at 164, and the resulting intermediate images are sub-
jected to up-sampling at 166 and 168 digital shear coordinate transformations at 168. Subsets
of those images are aggregated at 170, and selected resulting images are subjected to addi-
tional coordinate transformation at 172. This process continues until all of the projections
and intermediate images have been processed and aggregated at 174, to produce the image j.
The final coordinate transformation K_n//2 shown here at 1.72, only involves rearrangement
of pixels, when the sampling pattern for the pixel image / is Cartesian.
Fig. 14 is an algorithm for finding parameters of coordinate transformations for previ
ously described embodiments, using the properties of the Fourier properties of intermediate
images shown in Fig, 12.
Fig. 15 illustrates an algorithm for another embodiment of the present invention, over-
sampled ternary two-shear hierarchical backprojection. The embodiment of Fig. 15 is similar
to the embodiment of Fig. 10A, and the reference numerals from Fig. 10A are used where
appropriate. In addition to the steps of Fig. 10A, however, the embodiment of Fig. 15
includes a downsampling step 109 in the one before the last, which is the second level in
Fig. 15. Here too, the parameters of the various digital image coordinate transformation are
chosen to account for the view angles of the constituent projections of the intermediate im-
ages, and for the Fourier characteristics of the intermediate images, so that the aggregates of
the intermediate images may be accurately represented by sparse samples. However, certain

degrees of oversampling are used to improve the accuracy of subsequent processing, as will
be explained. If desired, for improved computational efficiency the downsampling step 109
can be combined with the second, x-coordinate shear comprising the two-shear digital image
coordinate transformation 108 (shown in Fig. 10B) producing a shear-scale transformation,
so that the processes 108 and 109 are together replaced by the process shown in Fig. 15(B).
Fig. 16 illustrates an algorithm for another embodiment of the present invention, over-
sampled ternary hierarchical shear-scale backprojection. The embodiment of Fig. 16 is
similar to the embodiment of Fig. 13, and reference numerals from Fig. 13 are used where
appropriate. However. Fig. 16 includes additional steps of upsampling 161 in me intermedi-
ate levels, and downsampling 169 in the level before the last level, in which the intermediate
images are upsampled and downsampled, respectively. Similarly to the embodiment of Fig.
15, in the embodiment of Fig. 16 when downsampling step 169 along the x coordinate.
ix U1,m follows an x-shear step 168 Sx',m, the two can be combined for computational effi-
ciency into a single digital image shear-scale.
Non-cartesian Sampling Schemes For Fast Hierarchical Backprojection Algorithms
The intermediate images in the different embodiments of hierarchical backprojection
illustrated in Figures 8,1.0, and 13 have a peculiar spectral support amenable to efficient non-
cartesian sampling. In particular, the underlying continuous domain image Itm in the Ith
level occupies a wedge in Fourier space, as seen for example in Figures 6, 7, 9, and 12. Multi-
dimensional sampling theory tells us that for images with a spectral support such as this,
sampling on a cartesian grid is less efficient than sampling on an appropriate non-cartesian
grid, Non-cartesian sampling can improve sampling efficiency by packing the copies of the
baseband spectrum more tightly in the Fourier plane. For an explanation of 2D sampling see
[?]. In particular, a quincunx sampling scheme reduces the sizes of intermediate images, and
therefore the computational costs of the algorithm, by a factor of almost 2. Digital image
coordinate transformation on periodic non-Cartesian sampling patterns can be executed ef-
ficiently using one dimensional shift-invariant filters. Likewise, all the methods previously
described for selection of the parameters of digital image transformations apply as well in the
case of such sampling patterns. Therefore, all the embodiments previously describe extend
to embodiments that use periodic non-Cartesian sampling patterns.
Variants of me embodiments of the present invention already described are applicable

to 3D backprojection of a variety of forms of 3D projections, including the X-ray transform
that arises for example in Positron Emission Tomography, and the 3D Radon transform,
which arises in magnetic resonance imaging.
The three-dimensional (3D) X-ray transform is a collection of integrals along sets of
parallel lines, at various orientations, in 3D. Each 3D X-ray projection is a two-dimensional
function that can be characterized by the 3D angle at which the lines are oriented. The block-
diagram for hierarchical backprojection for 3D X-ray data is similar to the ones previously
described, such as Figures 8. 10, 13, 15 or 16. The intermediate images in this case are
three-dimensional, not two-dimensional as in the previously described examples. Each in-
termediate image is sampled on a 3D sampling pattern that is sparsest in a direction that is
an average of the 3D angles of the constituent projections. As the algorithm progresses the
density of samples along this sparse-sampling (slow) direction increases to accommodate
the increasing bandwidth in that direction as explained by the Fourier analysis of 3D X-ray
projections. Consequently at every stage in the algorithm, before the images are aggregated,
each has to be upsampled along this slow direction. The extra dimension available in the 3D
embodiment also provides more coordinate transformations available for use in the various
steps in the algorithm, such as rotations in 3D. As in the 2D case, the parameters of these
digital image coordinate transformations can be chosen to account for the constituent view
angles and for the Fourier characteristics of the intermediate images. These digital image co-
ordinate transformations can be decomposed into a sequence of one-dimensional operations,
such as shears and shear-scales, as previously described. As in the 2D case, oversampling in
any subset of the 3 dimensions may be enforced to improve image quality.
A 3D radon transform projection is a one-dimensional function — a collection of in-
tegrals along sets of parallel 2D planes, parameterized by the displacement of the plane from
the origin. The view-angle of the projection is that of the 3D orientation angle of a vector
perpendicular to the set of planes. The block-diagram of the hierarchical backprojection of
3D radon projections is as in Figures 8, 10, 13, 15 or 16. In the first level the projections
are Radon backprojected onto the 3D image domain. These images are constant along two
dimensions, and therefore need be sampled only on the direction perpendicular to the two
constant directions. When groups of 3D radon projections are combined, the bandwidm of
the aggregate image may increase in one or two dimensions, depending on the view-angles of
these constituent projections. It is therefore necessary to upsample the intermediate image,
possibly in two dimensions, before coordinate transforming it and adding to other intermedi-

ate images. As in the 2D case, the coordinate transformations may be performed separably,
may be combined with the upsampling operation, and oversampling may be enforced on the
intermediate images.
In the various embodiments of the present invention, digital image coordinate trans-
formations and downsampling or upsampling operations may be performed by a sequence of
lower (one,) dimensional linear digital filtering operations. Furthermore, when the sampling
patterns used have some periodicity, these digital filters can be shift invariant, as will be de-
scribed in more detail. For computation efficiency, the digital filters can be implemented us-
ing recursive filter structures, or using an FFT, as is known. One way to determine preferred
values for the digital filters is using the theory of spline interpolation, as explained in M.
Unser, A. Aldroubi, and M. Eden, Fast B-spline transforms for continuous image represen-
tation and interpolation, IEEE Transactions on Pattern Analysis and Machine Intelligence,
Vol. 13, pp. 277-285, 1991; M. Unser, A. Aldroubi, and M. Eden, B-spline signal, process-
ing: Parr I-thf.ory, IEEE Transactions Signal Processing, Vol. 41, pp. 821-832, 1993; and
M- Unser, A. Aldroubi, and M. Eden, B-spline signal processing: Part II - efficient design
and application'!, IEEE Transactions Signal Processing, Vol. 41, pp. 834-848, 1993
The hierarchical backprojection methods of the present invention are applicable to
a wide range of tomographic imaging geometries, such as divergent beam geometries, in-
cluding fan-beam and cone-beam, with arbitrary source trajectories. The common fan-beam
geometry with a circular scanning trajectory is depicted in Fig. 17. The ray source moves
on a circular trajectory of radius D around an image of radius R. A fan-beam projection at
source angle B corresponds to line integral measurements along a set of rays parametrized
by fan angle y. TST is the distance along the ray of the source to the closest edge of the
image-disc and TEND is the distance of the source to the farthest edge of the disc.
Fig. 18 illustrates an algorithm for an embodiment of the present invention, hierarchi-
cal ternary rotation-based backprojection, applicable to fanbeam weighted backprojection.
The embodiment of Fig. 18 is similar to the embodiment of Fig. 8, and the reference numer-
als from Fig. 8 are used where appropriate. The embodiment of Fig, 18 differs from that
described in Fig. 8 in several respects. First, P = 4 x 3L projections, and the case of L — 2
is shown. Second, at the zeroth level, the initial intermediate images Id1,m are produced from
the fanbeam projections by weighted backprojection at 99, denoted here by W. Third, at
the last stage, four rather than two intermediate images are aggregated. Fourth, the sampling
patterns used in most of the early levels of the hierarchy are preferably chosen to be non-

Cartesian, as will be. explained. Here too, the last digital image coordinate transformations
only require reordering of image pixels, if a Cartesian sampling pattern is used for the in-
termediate images Id3m in the one-before-iast level Also, as in Fig. 5(B), the operations in
the zeroth and first level of the hierarchy can be combined, so that the intermediate images
Id2mi are produced by a direct weighted fanbeam backprojectkm of their three constituent
projections, or by other means.
In a selected number of levels, it is beneficial to modify the embodiment of Fig. 18,
by including additional weighting steps before and after the cascade of upsampling step and
digital image rotation (steps 106 and 10S in Figure 18), as shown in Fig. 19. The intermediate
image I2tm is weighted by spatially varying weights at 180 and 182, respectively before and
after steps 106 and 108. As will be explained, this weighting can be used to reduce the
sampling requirements of the intermediate images, and thus reduce the computation.
The sampling scheme of the intermediate images afects the performance of the algo-
rithm. The. desired sampling scheme is one that uses the fewest samples without losing image
information. Here too, the parameters of the various digital image coordinate transformations
and resampling operations can be chosen to account for the view angles of the constituent
projections of the intermediate images, and for the Fourier characteristics of the intermediate
images, so that the aggregates of the intermediate images may be accurately represented by
sparse samples, as will be explained. An alternative method for choosing these parameters
is based on intersection of particular sets of rays or curves, as will be described.
Figs. 20A-20C illustrate the progression of intermediate images through the levels
of the ternary hierarchical algorithm described in Fig. 18, for the case of source, angles
uniformly spaced at intervals Δp. The fans of the projections that make up intermediate
images in the algorithms are shown. At the zeroth level, each intermediate image Id1m is
made up of a single projection with fan oriented at β = 0, shown in Fig. 20(A). After the
first level of the recursion, each intermediate image is made up of three projections, with fans
as shown in Fig. 20(B) After the second level, each image is made up of nine backprojected
fans , as shown/in Fig. 20(C).
The intersection-based method for choosing coordinate transformations and sampling/
resampling patterns in the algorithm is illustrated in Figs. 21A amd 21B, for intermediate
image 73/m, with constituent fans as shown in Fig. 20(C). The sampling pattern for 73m
is made up of points that lie on the central constituent fan, which coincides with the one

shown in Fig. 20(A). As shown in Figs. 21(A) and 21(B), the sampling points for I3m are
determined by the intersection of half of the central fan with the extremal constituent fans on
the respective side of the central fan.
It is advantageous to modify the resulting sampling pattern by applying two con-
straints, to improve the accuracy of the back-projection, and reduce the computation require-
ments. First, the density of samples along the rays is limited not to exceed the desired
sampling density of the final image. Second, the sampling pattern is forced to contain on
each ray at least one sample on the outside of the image disk on each side. The plots in Fig.
22 displays the position of sampling points along a particular ray of the central fan at fan
angle y1. The sampling points lie at constant intervals in y', where Y is the fan angle of one
of the extremal fans shown in Figs. 2l(a)(b) and (c). The points that fall on the continuous
curve in Fig. 22 are original intersection points obtained using the method illustrated in Figs.
21A and 21B. The points on the broken curve are those modified using the above two con-
straints. ?Fig. 23(D)? shows an example of a sampling pattern obtained using this modified
intersection method. The fan shown is the central fan of the intermediate image for which
this sampling pattern has been produced.
As in the case of the previously described embodiment of the present invention, it
is advantageous to decompose the resampling and digital image rotations operations into a
sequence of one dimensional operations. The blocks marked f U in Fig. 18 represent the
upsampling of the intermediate image sampled on a central fan onto a finer set of sampling
points on the same fan. This can be achieved by separate upsampling along each ray of the
fan. For the cascade of upsampling and rotation marked by + U and K in Fig. 18, it is
conveniet to do the decomposition into ID operation jointly, as illustrated in Figs. 23(A)-
23(D). In each of the four panels, the two dashed circles represent the boundary of the image
of radius D and the source trajectory on the (larger)circle of radius R. The sampling points
are represented by small circles. The digital intermediate image with central fan (Fan,J and
the sampling pattern shown in Fig. 23(A) is to be resampled and rotated to the angle of the
central fan (Fan^) shown in Fig. 23(D), with the final sampling pattern shown in Fig. 23(d).
This is accomplished in the following two steps:
(i) Upsample the image, separately along each ray of Fana, to the sampling pattern shown
in 23(b), which is defined by the intersection of Fand with Fand. These points will
therefore lie on Fand, as shown in Fig. 23(c);

(ii) Upsample the image separately along each ray of Fand, to the sampling pattern shown
in Fig. 23(d).
These steps accomplish the combined operations of resampling and rotation, using 1D
upsampling operations.
The Fourier-analysis techniques described for the embodiments of the present inven-
tion in the case of intermediate images sampled on periodic sampling patterns are extended
to devise spatially varying sampling schemes in the divergent beam case. These techniques
are general enough to be applied to projections and back-projections on arbitrary trajectories.
in both two and three dimensions. For the case of backprojection on non-periodic systems of
lines or curves, the concept of local spectral support replaces that of spectral support. This
is illustrated in Fig. 24 for an intermediate image produced by the backprojection of a single
fan shown on the left side of Fig. 24. As will be explained, the local spectral of this contin-
uous intermediate image at the indicated point parametrized by r and 9 is the line segment
in the Fourier domain shown on the right in Fig. 24(b). The local spectral support at a point
of an intermediate image composed of the backprojection of multiple fanbeam projections
is shown in Fig. 25. On the left, the position of the point is indicated, as well as the range
of view angles of the constituent projections for the intermediate image. The bow-tie shaped
region on the right is the corresponding local spectral support.
This analysis of the local spectral support is used to determine local sampling require-
ments for the intermediate image. The resulting spatially nonuniform sampling patterns
are indicated by the dotted arcs in Figs. 26A and 26B for two typical intermediate fan-
backprojected images. The image boundary is indicated by the circle in broken line, and
only few points need be taken outside this boundary. This local Fourier sampling method
may be applied directly to find sampling schemes for arbitrary projection geometries over
lines, curves or planes over arbitrary dimensions.
The Fourier-based method for determining the sampling patterns for resampling and
coordinate transformation for hierarchical fan-beam backprojection can be further extended,
as will be explained, to sampling patterns lying on lines or curves other than central fans of
the intermediate images. Examples of such beneficial sets are illustrated in Figs. 27(A) and
27(B).

General Divergent-Beam Algorithms
The methods described for fanbeam backprojecticn, extend directly to other divergent-
beam geometries, including 3D cone-beam, which is one of the most important in modern
diagnostic radiology as will be described. Similarly to the fan-beam geometry, in which a
source of divergent rays moves on a trajectory (i.e., a circle) around the imaged object, in
the general divergent-beam geometry a source of divergent rays moves on a (not necessarily
circular) trajectory in 3D space around the imaged object. The detector surface measures
line integrals through the imaged object.
One embodiment of a hierarchical backprojection algorithm for a general divergent-
beam geometry can again be described by a block diagram similar to Figure 18, but mod-
ified in the following ways. First, at the zeroth level, the divergent-beam projections are
zero-backprojected at 99 with the appropriate divergent-beam single-view backprojection
W1 corresponding to a nominal "zero" source position, producing initial intermediate im-
ages. Second, because the trajectory of the source is rot necessarily circular, the constituent
divergent-beams of the intermediate images may not simply rotated as in the fan-beam ge-
ometry, but also translated, with respect to each other. The coordinate transformations K in
18 are selected accordingly. Third, depending on the presence of symmetries in the source
trajectory and positions, there may or may not be "free" coordinate transformations such as
the. pixel re-arrangement which replaces a n/2 rotation in the fan-beam algorithm.
As in the fan-beam case, the. initial intermediate images are processed hierarchically by
the algorithm. Analogously to the fan-beam case, intermediate images mat are close to each
other in position and orientation are aggregated, in order that the aggregated intermediate
image might be sampled sparsely. The 3D sampling patterns for the intermediate images
are determined by studying the structure of constituent weighted backprojected divergent-
beams that arc rotated and translated with respect to each other. One sampling pattern, as
illustrated in Figure 28, is where the samples are located along the rays of a divergent-
beam corresponding to the central constituent beam of the intermediate image. The sample-
spacing along each ray is chosen to ensure that all the constituent divergent-beams of that
intermediate image are sufficiently sampled along the ray. Alternatively, a more general
way is to use the local Fourier method to find the sampling pattern, as described previously
for the fan-beam case. Knowing how an intermediate image is composed of its constituent
projections, the local 3D Fourier structure of every intermediate image is determined. A 3D

local sampling matrix function at each point of the intermediate image is found that matches
the sampling requirements for the local Fourier support, as described in the fan-beam case.
This matrix function is then integrated (possibly numerically) over the image domain to
determine the position of the samples.
A separable method of up-sampling combined with rotation and translation onto a
new sampling beam is achieved similarly to the fan-beam case. It reduces the 3D coordinate
transformation and resampling operations to a sequence of ID resampling operations. As
shown in Figure 29 (a), each divergent-beam may be regarded as the intersection of a set o(
vertical planar fan-beams that are distributed in azimuthal angle, with a set of tilted planar
fan-beams at different elevation angles . The steps of the separable coordinate transformation
are as follows (as shown in Figure 29 (b))
1. The original divergent-beam is resampled onto the set of intersection of the rays of the
original with the vertical planes of the new divergent-beam . These points are therefore
located on the planes shared, by the final sampling points.
2. Steps 1 and 2 from the separable resampling in the fan-b&am case are performed for
each plane separately to resample onto the final set of points.
With a suitably efficient sampling scheme, the fast hierarchical backprojection algo-
rithm for divergent-beam can be expected to achieve large speedups with desirable accuracy.
As in the fan-beam case one might use a pseudo-beam that is modified (e.g., with the location
of the vertex moving farther away from the origin) in successive levels,
As will be evident to those skilled in the art, the methods of the present invention
are not limited to the examples of imaging geometries or specific embodiments described
herein. The methods are equally applicable to general problems of backprojection with other
geometries.
Figure 30 illustrates resampling-based backpiojection, as the sum of upsampled in-
termediate images formed by generalized backprojection of individual projections at source
positions /?p. It is similar to the rotation-based backprojection 4, but differs from it in two
respects. First, in the first step the pth (possibly processed) projection is subjected to a
weighted-backprojection 184 at the source position or orientation /?p, rather than at zero, as
is the case in 4. For example in the case of fanbeam projections, each projection is backpro-
jected ar the orientation of its source-angle. Second, before aggregation by addition at 188,

each initial intermediate image undergoes an upsampling operation at 186. This is neces-
sary, because the sampling pattern of each of these P single-projection intermediate images
is chosen to be an efficient and sparse pattern, so will usually be different for each projec-
tion, and often non-Cartesian. Before the intermediate images are aggregated, they need to
be resampled onto a common and denser grid. By itself, this structure offers no reduction
in operation count. However, it serves as a stepping stone to the introduction of some of the
fast hierarchical backprojection methods of the present invention.
Figure 31 is another embodiment of the present invention, a resampling-based hier-
archical backprojection for creating a pixel image f from a plurality of projections q1...qP.
Figure 31 is the binary hierarchical version of Figure 30. First, as in the non-hierarchical
case, each projection is backprojected 184 at its individual orientation, onto a sampling pat-
tern suited to that onentadon, producing an intermediate digital image. This is the zeroth
level of the hierarchy. In the first level of the hierarchy, these intermediate digital images are
upsampled at 186 to a denser sampling pattern common to selected pairs of images. This
upsampling will usually be to a non-Cartesian sampling pattern, but can be performed by a
sequence of one. dimensional resampling operations, as shown in Fig. 23, or in Fig. 29(b).
Selected resulting upsampled images are aggregated pairwise at 190, producing new inter-
mediate images. In the second level, the new intermediate images are again upsampled at
192, and aggregated at 194, producing new intermediate images. This process continues until
all intermediate image and projections have been processed, producing after the last aggre-
gation step 198 the final image /. As in the previously described embodiments, operations
can be combined within and across a level.
The sampling patterns at each stage in the hierarchy and the parameters of the up-
sampling operation in the embodiment shown in Figure 31 may be chosen by any of the
previously described methods. For example, in the case of fan-beam projections, one possi-
ble sampling pattern for a given intermediate image would lie on the points of intersections
of two fans: the first oriented at the central constituent source position; the second oriented at
an extremal constituent source position. Alternatively the sampling pattern and parameters
of the upsampling steps can be determined based on the Fourier or local Fourier proper-
ties and the view angle of projections included in the intermediate images. In particular, for
non-periodic sampling patterns, the local Fourier method described for the fanbeam rotation-
based algorithm can be used to find the sampling patterns: knowledge of how the projections
combine to form the intermediate images leads to the determination of the local spectral

support, which is used in turn to calculate the local sampling matrix function, which when
integrated over the image domain produces the sampling pattern for the intermediate image.
Fig 32 is the block-diagram for fast hierarchical reprojection. Reprojection is the pro-
cess of computing,to mo graphic projections from a given electronic image. The reprojection
algorithm is found by applying the process of flow graph transpositionto any block diagram
of a a backprojection algorithm, possibly with some change in weighing operations. In the
process, operations are replaced by their adjoint or dual. The block diagrams for reprojection
therefore appear similar to a reversed version of the corresponding one for backprojection.
The reprojection process described in Fig. 32 is one such embodiment of reprojection, ob-
tained from the backprojection algorithm described in Fig, 8.
in the first level , a copy of the input image f is preserved in the top branch of the
diagram as a top intermediate image, and in the bottom branch the image f is rotated at
200 by — n/2, producing a bottom intermediate image. In the second level, in the top half
of the diagram, the un-rotated to intermediate image is subject to three separate digital im-
age coordinate transformations at 202 some of which may leave the image unchanged, pro-
ducing three different top intermediate images. A similar process is applied in the bottom
branch, producing three bottom intermediate images. Each of the top and bottom intermedi-
ate images (six in all) then undergoes a process of decimation (low-pass filtering followed by
downsampling) at 204, producing new intermediate images. In the instance of the embodi-
ment illustrated in Fig 32 there are only 2 - 3L, with L — 2 view-angles, so the third level is
the final one in the recursive hierarchy. In the third and final level the intermediate images
are subject to separate coordinate transformations (2Q6)some of which may leave the image
unchanged, producing 18 intermediate images. The last step, which is not part of the recur-
sion is different: each intermediate image undergoes a reprojection at zero degree at 208.
Reprojection at zero degrees is equivalent to summing the vertical lines of pixels to produce
a one-dimensional, signal. These one-dimensional signals (518) are the output projections
produced by the algorithm.
The parameters of the digital image coordinate transformations in the algorithm are
chosen by the knowledge of the Fourier characteristics of the intermediate images. These
parameters are simply related to the parameters of the corresponding backprojection algo-
rithm. It is easy to see that since the reprojection block-diagram is a flow transposition of a
backprojection block-diagram, every branch of the reprojection block-diagram has a corre-
sponding branch in the backprojection block-diagram, and the coordinate transformations in

the corresponding branches are mathematical adjoints of each other. In the version of this
reprojection algorithm that corresponds to the two-shear backprojection algorithm in Fig 10,
the coordinate transformations in the second level (202) is an x-shear followed by a y-shear,
and the coordinate transformation in the last level (206) is an x-shear,fo!iowed by a y-shear,
and a fractional decimation in x (These three operations can be reduced to a shear-scale).
The parameters of these shears are the negative of the corresponding parameter used in the
backprojection algorithm. The parameter of the decimation in x is the same as that of the
upsampling in x in the first level of the backprojection algorithm. In the shear-seal^ version
of this algorithm, (corresponding to the shear-scale backprojection algorithm displayed in
Figure 13), the coordinate transformations in the second level (202) are shears in x (and the
parameter of each shear is the negative of tile corresponding- parameter in Figure 13). The
coordinate transformations in the final level (206) are shear-scales. The shear-scale used in
the reprojection algorithm is the mathematical adjoint of the corresponding shear-scale used
in the backprojection. The parameters of the decimation factors are also the same as the
corresponding upsampling factors in the backprojection.
Just like in the backprojection algorithms, oversampling of the intermediate images
in the algorithm can be enforced by first upsampling the images at the beginning of the
algorithm and downsampling them by the same factor at the end of the algorithm. Also.
these operations can be combined within or across a level.
Fig. 33 is the block-diagram of a decimation-based weighted reprojection. It it the
flow-graph transpose of Fig. 31 It shows the reprojection of the image onto a set of projec-
tions at 18 different source angles. Initially the given image is processed along two parallel
paths. In each path the image is subject to three parallel resamplings (210) onto three differ-
ent sparser sampling patterns. This resampling can Le performed in a separable way using
one-dimensional decimations (low-pass filtering followed by resampling onto a sparser grid).
The parameters of the filter are determined by the Fourier characteristics of the intermediate
image and the desired projections. Local Fourier analysis of the desired projections is used
in the case when projections do not line on parallel straight lines. In the final level of the al-
gorithm each intermediate image is again subjected to three parallel resamplings (212) onto
sparser sampling patterns. Finally a weighted projection is performed on the image to pro-
duce the projection p0. This involves a weighted sum of the pixels of the image to produce a
one-dimensional projection.

Implementations and Experimental Results
Preferred embodiments of the present invention were implemented in C programming
language and tested in numerical experiments on a Sun Ultra 5 workstation with 384 MB
RAM. The test image was the Shepp-Logan head phantom(a standard test image used in
numerical evaluations of tomographic algorithms) By varying the parameters of the algo-
rithms a tradeoff can be made between accuracy and computational cost. Accuracy refers to
the quality of the reconstructed image, Though visual quality is not easily quantifiable, we
measure the error between the reconstructed image and the original from which the radon
transform was numerically computed. The measure of error used is the relative root-mean -
square. error (RRMSE). The RRMSE in reconstructing an N x N image /[m2, rnL] from the
tomographic projections of f[rn,2,m1] is calculated as follows:

For parallel-beam data, the test image was of size 256 x 256, the number of view an-
gles was 486, and the Shepp-Logan filter (the ideal ramp filter multiplied by a sine function)
was used to filter individual projections. Fig. 30 displays the RRMSE error versus the run
times for the two algorithms at various values of the oversampling parameter 7 between 0 75
and 1.0. The two-shear algorithm is represented by the circles and the shear-scale algorithm
by the squares. Each algorithm is run using two types of filters — a third-order (dashed line)
and fifth-order (solid line) spline filter called MOMS 16. For each flavor of tire algorithm,
as 7 is decreased the error of the algorithm decreases and the run time increases. The plot
points that are not connected to any other are the non-oversampled versions of the algorithms
represented by the connected points. In comparison the run time of the conventional algo-
rithm, using linear interpolation, is 14s and the RRMSE error of its output image is 0.04S6
(worse than the fast algorithms displayed here).
Some sample images from the output of the algorithms, for parallel-beam data, are
displayed in Fig. 35. Columnwise from left to right, they are output images from the con-
ventional backprojection, the two-shear and the shear-scale algorithms. An oversampling of
7 = 0.82 was applied to the two fast algorithms. The lower row of images displays in detail
a section of the corresponding images in the upper row.

The fast hierarchical algorithm for fan-beam geometry was successfully tested on the
512 x 512 2D shepp-logan phantom. The acquisition geometry considered was with a source
radius D = 1.06 x N = 544, 972 source angles, and 1025 equiangular detectors. The
regular and oversampled version of the fast algorithm was implemented using a variety of
interpolation methods. The resulting reconstructed images were compared to a conventional
fan-beam algorithm that used linear interpolation. In all the experiments the projections were,
filtered with the Shepp-Loga.n filter.
Figs. 36(A) and 36(B) display the reconstructed images from the conventional in Fig.
36(A) and the fast algorithms in Fig. 36(B). The result of the fast algorithm is comparable
to that of the conventional algorithm. The point spread functions of the fast algorithms are
comparable to that of the conventional one. Fig. 37(B) displays the PSF of the conventional
algorithm, Fig. 37(C) displays the PSF of the fast algorithm with linear interpolation and
7 = 0.4 oversampling. Figure 37 (c) compares slices through the x-axis of the psfs of the
conventional algorithm, the non-oversampled and the 7 — 0.4 oversampled fast algorithms.
The similarity of the PSFs confirms the comparable image quality in the fast algorithms.
The sampling scheme used was the fan intersection-based method, without the en-
hancements suggested later. In addition to the shortcomings of this scheme mentioned pre-
viously, for reasons of simplicity of implementation, numerous sample points nor needed for
the correct operation of the algorithm were used in the embodiment tested in the experiments.
Despite these inefficiencies, for N = 512 and D = 1.41 N, the ratio of (data-dependant) mul-
tiply operations in the conventional algorithm to the fast (linear) one was 6.4 and the ratio of
addition operations was 3.0. Note that the geometry used in this experiment, with a source
very close to the origin (D = 2.06 R) is particularly challenging for this un-enhanced imple-
mentation, because of the high density of samples near the source vertex. In most practical
systems, the source is further away, reducing these effects. Furthermore using the alternative
sampling schemes discussed earlier, will lead to much higher speedup factors.
Detailed Description of Backprojection
and Reprojection Algorithms Used In The Present Invention
Backprojection is the process used to create images from processed projections. Re-
projection is the reverse process, used to compute projections from a given image. Both
operations are used in image reconstruction from projections, as seen in Fig. 1. Conven-

lional backprojection will be described first, followed by a description of the backprojection
and reprojection methods of the present invention. The description will be given first for the
case of parallel-beam projections of a 2D image, and then for the more general 2D and 3D
geometries,
The classic and preferred method used to estimate an image from its projections is the
filtered backprojection (FBP) algorithm. The FBP involves first filtering the projections with
a so-called ramp or modified ramp filter and then backprojecting those filtered projections.
Additional pre-processing steps may also be applied to the projections before backprojection.
A processed projection taken at view angle 0P, wherep = 1, ....P is the projection index, will
be denoted alternately by q(t, p), or qP or q0„, depending on whether the dependence on the
variable t or 8 needs to be shown explicitly. For brevity, processed projections will be called
projections hereinafter.
The FBP reconstruction f from a set of P projections (at the angles specified by the
components of the length-P vector 0), is described by:

where
Definition 0.1 the backprojection operator is defined by

The 0(N2logP) coordinate-transformation-based hierarchical backprojection algo-
rithms of the present invention are based on a hierarchical decomposition of backprojection
in terms of coordinate transformations. The details of several preferred embodiments are. de-
scribed with different choices of coordinate transformation, concluding with the most general
coordinate transformations.
In particular, the 0(N2 log P) rotation-basec hierarchical backprojection algorithms
of the present invention are based on decomposition of backprojection in terms of the rotation
of images.
The rotation-by-0 operator K.{0) which maps an image f(f) to its rotated version


As seen in Fig. 4, Lhe backprojection in equation (3) may be rewritten as fallows. Here
qv refers to Che projection whose view-angle is 0p,

One embodiment of hierarchical backprojection stems from the fact that the cumula-
tive result of several successive rotations is still a rotation. In particular,

It follows from Equation (6) that with P = 2L for some integer L, the block diagram in Fig.
4 can be rearranged into a hierarchical tree structure as shown in Fig. 5. The intermediate
image in the mth branch of the Ith level is denoted as I1m. In the initial level

intermediate images are produced by backprojecting individual projections, and in subse-
quent levels, i.e., for / — 1, 2, 3,..., log2 P


the intermediate images undergo coordinate transformations (rotations, in this embodiment)
and are aggregated by addition. Fonany set of projection angles 0i, i— 1,..., P, the interme-
diate rotation angles δ1,m can be chosen to guarantee that the structures in Fig. 4 and 5 are
equivalent, so mat the hierarchical backprojection algorithm depicted in Fig. 5 produces the
desired image f.
In fact, because there are many more intermediate rotation angles (free parameters)
than view angles (constraints), there are many degrees of freedom in the choice of the δtm.
Also, for the digital images used in practice, digital image coordinate transformations are
used, whose accuracy and computational cost depend on the choice of sampling patterns and
implementation. These various degrees of freedom, or parameters, are used in the present
invention to reduce the computational requirements of the hierarchical backprojection algo-
rithm, as will be described.
Definition 0.5 The set NiiW, of indices of constituent projections for- an intermediate image
IltTn is the set of indices of projections for which there is a path from the input in Fig. 5 to the
intermediate image Iim.
For example, as seen in Fig. 5, It is easily shown, using
equation (8) or upon examination of the block diagram, that

Thus NiiTn lists the projections that make' up intermediate image ll>m. If, instead of being
processed by the method described in Fig. 5, these same projections indexed by the set N[iTn
were directly backprojected together at their respective view-angles, this would produce the
virtual intermediate inwge I^m:

Upon examination of the block diagram in Fig. 5, it can be seen that the relative angle
between projections in an intermediate image is preserved in the final reconstructed image
/. In other words, for all intermediate images in the algorithm, i.e. VZ, m


The intermediate rotation angles 6irm of the hierarchical algorithm are completely de-
termined by Q7>m as follows- It follows easily from the definition of Ntirn or Equation (9) chat
iVi+1>rn = Nt,2m~i U M,2m- and consequently by the Definition /.

Equations (S). (10) and (11) imply that

The intermediate rotation angles ( are chosen to reduce the bandwidth of the intermediate images. Images of small bandwidth
can be represented by few samples, and consequently reduce the computational cost of the.
whole algorithm. The bandwidth of these intermediate images is determined by understand-
ing the algorithm in the Fourier domain and tracking the evolution of the spectral support of
the intermediate images through the hierarchy.
Tomographic projections in the Fourier domain
The parameters of the digital image coordinate transformations are chosen to account
. for the view angles 0P of the selected projections, as has been described. The parameters
are also chosen for their effect on'the Fourier characteristics of the intermediate images.
These considerations are also used to determine which intermediate images are selected for
aggregation together.
Key to interpreting the backprojection algorithm in the Fourier domain is the projection -
slice theorem that relates the one-dimensional Fourier transform (J^i) of a projection Vgf
with the two-dimensional Fourier transform (T2) of the image /. The projection-slice theo-
rem 4 says that
The backprojection algorithm of Fig, 4 can be interpreted in the Fourier domain. The
backprojection-at-zero operator produces an image whose spectral support is limited to the

horizontal frequency axis u)x:

It is well known that rotating an image in the space domain results in the rotation of
its Fourier transform by the same angle

It follows that the function of the baekprojection operator in Equation (5) is to rotate
the spectral component of each projection to the appropriate angle in the reconstructed image
and add them all together.
Assuming,projections with one-sided bandwidth (in the t variable) equal to IT, the
typical virtual continuous intermediate image 7;;Tn in the hierarchical algorithm therefore has
a spectral support of a wedge in a range of angles determined by the constituent projections,
as shown in Fig. 6. In particular, let I^m be an intermediate image in the block diagram
of the algorithm shown in Fig. 5(A), with Ni>m = {b,b 4- l,...,e}7 i.e. the view-angles
of the constituent projections of this image are {6V : p = t, b + 1,..., e}. Because of the
relationship in Equation (10), the spectral support of IiiPt is just a rotated version of that of
Ii Fig. 6 is [7 — /?, 7 + j3] — [9b — oniTn, Bt ~ a(tm]. The bandwidths of Ii.m in the first and second
coordinate, fl3 and O/, respectively, are indicated on the bounding rectangle for the wedge,
which is shown by broken lines in Fig. 6. Thus, the choice of or!]T71 provides a way to control
the bandwidths of Ii>m.
Sampling theory dictates how such a continuous intermediate image IiiVl may be repre-
sented by a samples on a Cartesian sampling pattern. In particular, a continuous image with
bandwidths of Qs and Q/ in the first and second coordinate may be sampled on a Cartesian
pattern aligned with the coordinate axes, with a sample spacing of A, and A/ in the first and
second coordinates such that Afi of samples required to represent the continuous im£.ge, the quantity 1/(A5A/) > n2QsQf,
which is proportional to the area of the bounding rectangle, should be minimized. The rota-
tion angle that minimizes mis area, and therefore minimizes the sampling requirements for

intermediate image i}im, is
Selecting Parameters for Coordinate Transformations
Equation (15) is used in Equation (12) to determine the optimum rotation parameters
8im for the digital image coordinate rotations in the rotation-based hierarchical backprojec-
tion embodiment illustrated in Fig, 5. The indices b and e are determined as corresponding
to the smallest and largest, respectively, in the set NiiTn. With this choice, the intermediate
image ILw = JC(alm)Iiim, and its bandwidths in the first coordinate (slow bandwidth) and
second coordinate (fast bandwidth) are, respectively,

To further minimize the bandwidth, the differences 6&—8h should be minimized, which
is achieved by selecting to aggregate at each step intermediate images produced from pro-
jections with the least maximum angular separation between their view angles.
In addition to rotation angles and selection of which intermediate images to aggre-
gate, it is necessary to choose appropriate sampling patterns for the digital image coordinate
transformations of the various intermediate images in the algorithm. These are chosen us-
ing the sampling requirements, which are in turn determined using the spectral supports and
bandwidths of the intermediate images.
Jn particular, for initial intermediate images formed by the backprojection of a single
projection, Ilt,n — BoqSm as in Fig. 5, the slow bandwidth Qa(Ii,m) = 0, and a single sample
along the x2 coordinate suffices. The slow bandwidth will be larger and more samples along
the x-i coordmate will be needed for initial intermediate images formed, as in Fig 3-b, by
backprojcciing more than one projection.
The sum of two virtual intermediate images (lijin = Ii-i2m-\ + h-x^m) an^ their cor-
responding intermediate images are illustrated, in the space and frequency domains, in Figs.
7(A)-7(C). Figs.' 7(A) and 7(B) show the Fourier-domain support of / and /. TFig. 7(C)
shows the space-domain sampling scheme of IitmJi~itzm aud J/_li2m_i. The sampling points
are at the intersections of the horizontal and vertical lines in the space domain. The aggre-

gated intermediate image 7!im, having a greater slow bandwidth than each of its components,
requires a denser sampling pattern. Conversely, intermediate images at earlier levels of the
hierarchical algorithm require a sparser sampling pattern, leading to computational savings.
Computation Cost and Savings
The computational cost is readily estimated. Assuming equally P projection view-
angles in [0, IT), with an angular spacing of A^ = ir/P, and NitTn as specified in (9), equation
(16) (using b = 2'-l{rn - 1) + 1 and e = 2l~1(m - 1) + 2l~l) implies that 0,(7^) =
7rsin(7r(2'-1 - 1)/(2P)) and ,Q/(/(>m) = TX. The size (number of points.) of the digital
image lfm in the Ith level is therefore

Given that the cost of rotating a digital image of S pixels is 0{S) and size(/,rfm) =
0(N-2'/P). the complexity of the algorithm is

For the typical P = O(N), this is 0(N2 log N), which is much more favorable scaling
of the cost with the size of the image, than the 0(N3) of conventional backprojection.
Improved ternary rotation-based hierarchical backprojection
It is easy to see that the rotation of a sampled image by certain special angles —
0,±TT/2 and ix — are computationally inexpensive operations because they involve no inter-
polation but, at most, a mere rearranging of existing pixels. The computational efficiency
of the algorithm of the present invention can be improved by incorporating these free opera-
tions into it. The binary algorithm (Fig. 5) can be modified such that three, not two, images
are added ar every stage. The center image of each such triplet is rotated by 0 radians — a

free operation. To use the free rotation by — 7r/2, it is included in the final stage: one of the.
constituent images is rotated by — 7r/2 and added to the other tin-rotated one. This partic-
ular combination of binary and ternary stages results in a set of view-angles of size 2 x 'iL
for some integer L. Though the algorithm can be tailored to arbitrary sets and numbers of
projections, the description and analysis of the algorithm can be simplified by assuming that
they number exactly 2 x 3L and are uniformly distributed as follows:

Hence, the view-angles can be divided into two sets of 3L each, one set centered around
9 = 0 and the other centered around 6 — n/2. The block diagram of this ternary algorithm is
shown inFig. 8 for the particular case of L — 2, i.e., for a set of P = 2x32 = 18 projections.
The digital intermediate images lfm are the sampled versions of the underlying continuous
intermediate images J;|m, In the block diagram, blocks marked Z30 represent the zero-angle
backprojection operator. Blocks marked Kg represent the digital image rotation operator. As
will be explained, the sampling periods and rotation angles for the digital images and for the
digital image rotations are chosen to minimize computational cost.
Selected aggregate intermediate digital images can be stretched to produce new in-
termediate images. The blocks labelled -\y Uiiin in Fig. 8 represent up-sampling along
the second (slow) coordinate, which are necessary to accommodate the increasing band-
width of the intermediate images as the algorithm progresses. For the configuration of
view-angles in (18) and (22), the bandwidth of the intermediate images in the Ith level are
^.i(Zi,m) — Trsin. (A(3/_1 - l)/5s) and 0/(Ijjm) = n. The normalized spectral supports at
key points of the algorithm, numbered (1),...,(9), are displayed in Fig. 9. As the algorithm
progresses (as I increases), the slow bandwidth increases, and the required slow-sampling
period, A5 the various upsampling blocks can be chosen to satisfy these requirements.
Section 0.0.0.1 describes how the separable rotation of discrete images involving two-
shears is incorporated into the backprojection algorithm. In order to improve the quality
of the backprojected image, a systematic oversampling of the intermediate images is intro-
duced. Section 0.0.0,1 describes the algorithms modified to include this oversampling.

Intermediate rotation angles for coordinate transformation can be chosen as follows.
For general L, the hierarchy consists of (L+l) levels. In the zeroth level, each projection qg
is backprojected-at-zero (B0) to produce a corresponding image B0qg. The images are then
grouped into threes and combined co produce a third as many images. The groupings at level
/ are defined by the following relationships. In subsequent levels, / = l,,..,L,

For the configuration of view angles in (18), the optimal intermediate rotation angles are
£(,3m-i = 0 for Z — 1, 2, ...,L , ££,+1,1 = 0, and SL+Xt2 — —TT/2. Upon examination of Fig.
10, it is easy to see that the indices JV;jm of the projections that constitute image Il follows:

For the final stage we make use of the free rotations by ir/2 and 0 radians.
0.0.0,1 Coordinate transformation can be based on shearing, as follows. The digital
image rotation of the intermediate images can be replaced by a sequence of two digital image
shears as shown in Fig. 10(B), where shears in the x and y coordinate are defined as follows

Definition 0.7 The x-shear andy-shear operators are defined by
(Sx(a)f)(x) = f(S:x)
iSv{a)f){x) = f(S«x)
where the x-shear andy-shear matrices are S° - [\ f ] and S^ - [£ §], respectively
Tliis alternative embodiment is derived from the two-shear factorization of the rotation
matrix: KB = ^MflST;sin9c09X where Se = [co0s9 i/^*]. .This decomposition implies
that in the case of digital images, using perfect bandlimited interpolation and a sufficiently
bandiimited digital image as the input to the two-shear digital image coordinate transforma-
tion, the twice-sheared image is simply the image rotated by 8, and effectively downsampled
in x and upsampled in y by a factor of 1/ cos 6. To correct these changes in sampling pat-
tern would require non-integer resampling every time a two-shear operation is performed.
Instead, fay sampling the intermediate images differently, only the cumulative effects need to
be addressed for the final image.
A convenient choice is to correct the resampling once at the beginning of the hier-
archy, by upsampling the initial digital intermediate images in x and downsampling in y.
When the initial intermediate images are produced by single-view backprojection as in the
embodiments shown in Fig. 10, the downsampling in y is free, because it only involves
* initial backprojection onto a sampling pattern with a different y density. It is therefore not
shown explicitly in the block diagram in Fig. 10. Upsampling in x in the first stage avoids
aliasing problems in later stages. For each of the initial images, the cumulative effective
downsampling/upsampling Factor l/rL=iCOB^ to tne ^^ image is calculated across all
two-shear transformations on the path co the final image, and the reverse operation is applied
at the initial image as described. It follows that each of the initial images If ,p — 1,. ., P is
resampled with different x and y sampling densities.
So the resulting algorithm, as shown in Fig. 10(A) involves first resampling all the
images B0qBp,p = 1,..., P by the appropriate factors (depending on 8P) and then performing
two-shear digital image coordinate transformations as shown in Fig. 10(B). The two shears
are performed taking into consideration the changing sampling pattern of each intermediate
image.

Shear-scale algorithms
The shear-scale-based hierarchical backprojection algorithm is based on the definition
of backprojection as the sum of single-projection images that are scaled (stretched/contracted)
and sheared.
Definition 0,8 The x-shear-scale operator C(a) is defined by
where the x-shear-scale matrix
So, Equation (3) may be rewritten as

The following are some results describing the effect of accumulated shears and shear-
scales that can be easily shown.

Similarly to Equation (24), Equation (25) can be inductively applied to prove that a
sequence of n shear-scales is also a shear-scale. Using this property, Equation (23) can be
rewritten in a hierarchical algorithmic structure. Fig. 11(A) is the block-diagram representa-
tion of Equation(23) (for P=4), and Fig. 11(B) is its hierarchical equivalent.
The requirement of equivalence of the two structures sets up a system of equations
between the given projection-angles 8P and the unknown shear-scale parameters 5^, whose
solution always exists.
It is computationally beneficial (in terms of operations or storage space) to use an
integer scale-factor rather than a non-integer one. For example, it is beneficial to use x-





h+l,m — 'SI(tt|i3m-2)i'[|3m-2 + //,3m-1 + 4(ft!,3m)^,3m (-&)

initial level of the hierarchy depend on the set {ai]Tn} as follows :

where /i(p, 1) = [p/3'_1] describes the path in the hierarchy from the pth projection in level
1 to the root.
By examination of the upper half p the sequence of shears backwards from the last level, it can be seen that in this top half of
the hierarchy the underlying continuous intermediate image Il virtual intermediate image Ii
Similarly to the derivation of Equation (12), it can be shown that the intermediate
shear-factors depend on the {PitTn} as follows :

The freedom in the choice of the shear coefficients ai of the 0 parameters, which can be used to minimize the bandwidth of the intermediate pro-
jections. This will reduce their sampling requirements, and improve the computational effi-
ciency of the algorithm.
Let the spectral support of the intermediate image Ii^m-x be denoted by Wijm-%, and
its y bandwidth (i.e., slow bandwidth) be denoted by D,y(Wii3m-.i). The optimal Pi^m^
minimizing ^(Wi^-i) is then, for i = 0,1, 2

If the view angles are uniformly distributed, it is found that /3'3m_1 ~ PUim- An
advantageous choice is then /3,*3m_1 — A*+i,m> yielding a.\-im_l = 0, which eliminates about

one third of the required shear operations.
The set {/5;iTn} is determined by tracking the spectral support of the intermediate im-
ages in the algorithm, backwards from the final to first level of the hierarchy. Following from
the slice projection theorem, it is clear that the spectral support of I{m is

For the sake of sampling requirements, the spectral support W/.,,, is characterized by its end
points denoted as Ei,m.,

and shown in Fig. 12. The upper and lower band-edges of this set EisTn are denoted by
Efm and E{ respectively. So if Ni;m = {6, b + l,...,e}, I? Efm ~ (cos 9b, sin 6b). Because of the ternary combination rule, the upper, middle and lower
third sets of these spectral-support .end-points are denoted as Ef+l m,E}+1 , and Ef_l:m i.e..
Ell+1 m = Ei:im-i. And following from (30) and Fourier theory about affine transformations.
EitW — Sv(i3i_ni)Eitm. The optimal shear factors are then

The algorithm FlNDALPHAS shown in Fig. 14 finds the optimal shear-factors while traversing through the hierarchy down from level L to level 2 in the shear-scale back-
projection algorithm (Fig. 13). The inputs to the algorithm are the sets of end-points of the
spectral support of the images at the start of the (L + l)-th level: £?£+].,m for m — 1, 2. In
particular, in the case of uniformly distributed angles, EL+I,I = #L+I,2 = {(cos dpy sin 6p) :
p=l,2,...,P/2}.
Finding Upsampling Factors
The y-upsampling factors {UiyW} in the algorithms and different embodiments of this
invention determine the sampling density of the intermediate images in the slow direction.
These upsampling factors can be chosen to meet the sampling requirements of the interme-

diate images, while being restricted to integer values. This restriction to integers has com-
putational advantages, and may be expanded to that of rational numbers with low-valued
denominators for the same reason.
The problem of choosing these factors is slightly simpler in the case of the shear
scale algorithm than the two-shear case. In the former case, examining Fig. 13 it is easy
to see that the slow sampling interval of the intermediate image J/>m is Ylfi-t £Ax/,',/im)-
where ^(/', I, m) = fm/3'"'] describes the path from 7;,m to the final node of the algorithm.
Sampling theory requires that slow direction sampling interval be as follows :

where VLy{Wiim) is the y-bandwidth of I[>m. The computational cost of the whole algorithm,
given a set of upsampling factors hi == {UiiTn eZ:l = 2,3,.--, L;m = 1, 2, ..., 6 ■ '3L~1},
can be. shown to be
where the constants c(m denote the relative computational cost of digital coordinate trans
formation applied to intermediate digital image If These constants are determined by the
order of digital filters used, and by the particular implementation of the coordinate trans-
formation, as will be described in the discussion of efficient implementations of coordinate
transformations. The best set U to mixiirruze J(]U) subject to the constraint Equation (35),
can be. solved by dynamic prograraming or a comprehensive search. For a given set of view
angles, the best set of upsampling factors can be precomputed and stored in a lookup table.
In the case of the two-shear rotation algorithm, the problem is complicated slighdy by
the fact that the use of two-shear rotation causes a defacto fractional up and down sampling
of intermediate images. In particular, the slow-sampling periods A/>m in adjacent levels of
the algorithm are related as follows; As+ ,|m'3' = A',m/(C/j,m/c/)T7l) where


So the constraint on the slow direction sampling interval becomes as follows:

where n;/(M/!m) is the y-bandwidth of I[im. The computational cost given by Equation (36)
now can be minimized subject to the constraint in Equation (37).
Oversampled Versions of the Algorithms
Oversampling can be incorporated into this algorithm to improve the accuracy of the.
reconstruction. One way to do so is uniformly over all intermediate images in the algorimm:
whenever an interpolation is performed in the algorithm, the image being operated upon is at
least oversampled by some predetermined constant 7. In particular, the ratio of the Nyquist
rate to the sampling frequency (in both slow and fast directions) of every intermediate image
that is subject co a digital coordinate transformation or resampling should preferably be less
than 7. This oversampling is preferably incorporated while modifying the algorithms of the
present invention as little as possible.'
Oversampling in the slow direction
In the previously described algorithms of the present invention, the sampling fre-
quency in the slow direction is controlled by the upsampling factors Ui>m. This proves to
be a useful tool for maintaining the oversampling condition in the oversampled versions
of the algorithms, and therefore results in no alterations to the structures of the algorithms
of the. present invention, The upsampling by a factor of 1/7 is achieved by simply modi-
fying the constraint on the slow-direction sampling interval by the factor 7, i.e., requiring

Oversampling in the fast direction
In the fast direction, the computationally inexpensive integer upsampling is not used to
control the oversampling, so the algorithm is modified to involve fractional upsampling. In
the two-shear hierarchical backprojection algorithm, such a fractional upsampling by factors

{Uitm '- rn. = 1,2,..., 2.3L} in the slow direction is already included in level I = 1. We there-
fore simply increase the upsampling factors UilTn to incorporate this oversampling and then
downsample the image after the last digital coordinate transformation has been performed,
to return the image to the desired sampling scheme (where A/ = A,, = 1.0). This modi-
fication co the two-shear algorithm therefore involves only one additional level of fractional
resampling. This x coordinate resampling is combined with the digital 2-shear in the Lth
level for improved computational efficiency. The block diagram of this algorithm is in Fig.
15.
In the shear-scale hierarchical backprojection algorithm, the x-upsampling operations
are combined with the ^-shear-scales (C(a) at the beginning of the algorithm) to avoid an
extra level of resampling. It is easy to see that :c-upsampling is simply a special case of
x-shear-scaling, with the shear factor set to 0. The combination of these two operations is
still an x-shear-scale: C[a\, cr2)C(l/t/, 0) = C{a\/U, cr2/U). The downsampling at: the end
of the algorithm is combined with the s-shears in Lth level ( effectively shear-scales of the digital image. It can be shown that an x-shear by a followed
by a z-downsampling by U is effectively an x-shear-scale: C(U, 0)Sx(a) ~ C(U, 1 -I- Ua).
The exact values of these upsampling and downsampling factors is determined both by
the parameter 7 and the spectral structure .of the intermediate images. In the fast direction,
the oversampling condition is that the fast direction (y coordinate) sampling interval satisfy
Aj'm
Consequently, the upsampling factors Ui fied from U^m, those of the non-oversampled two-shear algorithm, as follows:

In other words, the upsampling factors in the first level are modified to ensure that
all the intermediate images in the algorithm are oversampled according to the parameter
7. Consequently, after the last rotation, the images have a fast-sampling interval of 1/r/;
and, therefore, at the Lth level, they have to be downsampled by 77 to return them to a unit-
sampling scheme.
4

It muse be noted here that because the input projections are assumed to be sampled ex-
actly at the Nyquist rate, the oversarnpling condition in the fast direction will not be satisfied
for the first-level images, B0qsp,p = 1,.... P.
The block diagram of this ternary oversampled shear-scale-based algorithm is in Fig.
16.
Collapsing sequences of similar operations
Whenever there is a sequence of two operations operating on the same single coordi-
nate, they may be combined for improved computational cost and resampling accuracy. The
following is another example, in addition to the previously mentioned ones of combining x-
shears and x-upsampling/downsampling, or x-shear-scales with x-upsampling/downsampling.
The non-oversampled two-shear algorithm of Fig. 10 involves a y-sheai followed by
an a'-shear of four out of six of the intermediate images in the Lth stage. We incorporate the
final downsampling of the image in the oversampled version of Fig. 15 by collapsing the
sequence of two operations — the ^-downsampling followed by the x-shear of this stage —
into a single one. This leaves the length of the cascade of interpolations unchanged from the
non-oversampling case.
Hierarchies of Arbitrary Radixes/Branching factors
All these algorithms can be easily modified for the case when the set of view-angles
is not of the form 2 + 3L. Though the preferred embodiment is for all the branches of the
hierarchy to involve the aggregation of triplets of intermediate images, or use rotations by
±7r/2 or 0, arbitrary numbers of intermediate images may be combined at each stage.
Given an arbitrary number (say M) of intermediate images at a particular level of the
algorithm, 3 x [M/3\ of mem may be combined in groups of three (where [x\ is the largest
integer less than or equal to x). If the remaining number of intermediate images is two then
they may be aggregated as a pair to produce an image at the next level. If there is only
a single intermediate image remaining it may be passed on, without alteration, to the next
level of the hierarchy.

The branching factor of the hierarchy (the number of branches that aggregate at a node
of the hierarchy) may be altered, not just to accommodate arbitrary numbers of view-angles,
but also to reduce the depth of the hierarchy and thereby improve image quality. In that case
it may be useful to aggregate images not in pairs or triplets but in larger groups.
The previous prescriptions for the parameters of the coordinate transformations may
be easily extended to nodes with arbitrary branching factors. For example in the. rotation-
based algorithm where /;+i,m = J2t=i £(^,m;)^!,m,-, the relations described in Equations
(10) and (15) still hold, and the rotation angles therefore are prescribed by ai^rm- In the case of the shear-scale based algorithm where //+lr7Tl — X^fi S-x{ai,rn.i)Ii,mi.
Equation (29) still holds. Extending the notation E\tTn to refer to the ith set of the M sets of
projections being aggregated, one obtains an equation for a* identical to the right-hand-
side of Equation (34).
0.1 Hierarchical Algorithms Based On Other Image Transformations
The back-projection equation (3) may be written in the form

for any matrix Ag, as long as the first row of Ag is [cosesme]. The freedom to choose
arbitrary values for the remaining two entries of Ag allows for flexibility in the design ,pf the
coordinate transformations used in the hierarchical algorithm. Matrix Ae can be factored as
Ag — J|;t=i A(6i) for some parameter vectors 6i that are related to 0, but are not completely
determined, owing to the freedom in the two bottom entries of Ag. This factorization can
be used to derive a hierarchical decomposition of the backprojection equation (40) with
a corresponding block diagram such as Fig. 5, with the coordinate transformation steps
denoted by K.{8^m) representing the image coordinate transformation operators defined as
(1C(6i,m)f)(z) = f(A(5t,m)ti).
Clearly, the specific embodiments of this inversion described herein are special cases
of this more general choice of matrix A(g and its factorization, and the associate coordinate
transformation. The effect of these coordinate transformations on the Fourier spectrum of
the intermediate images is analyzed similarly to the cases already described, because the

effect of an affine transformation by matrix A in the spatial domain is an affine transfor-
mation by matrix A~T, i.e., the inverse-transpose of A, in the frequency domain. Similar
considerations can therefore be used to select free parameters in the transformations with the
goal of reducing the computational requirements. Thus, the class of digital image coordinate
transformations used in the hierarchical backprejection algorithms of the present invention
includes many other transformations in addition to those described for specific embodiments.
Furthermore, because matrices A[5^m) can be factorized into a product of triangular matri-
ces, the coordinate transformations can be performed as a cascade of single-axis coordinate,
transformations, if desired.
0.2 Efficient and accurate implementation of Digital Image Coordinate
Transformations and Resampling
The accuracy and speed of the hierarchical backprojection and reprojection algorithms
of the present invention depend on the specifics of the implementation of the various digital
image coordinate transformations and resampling operations. Improved accuracy requires
in general high order filters or interpolations, which usually increases the computation cost.
The cost of high order filtering or interpolation can be reduced by decomposing the op-
erations into lower order operations. Additional reduction in computation and/or memory
requirements is obtained if the filters used are shift invariant. High order finite impulse, re-
sponse filters can be implemented efficiently using low-order recursive filter, or by using the
fast Fourier transform (FFT):
In particular, if a separable representation basis is assumed for the continuous image,
a digital y-shear of the digital image can be achieved by fractional delays of the individual
vertical lines in a digital image array. Likewise, a digital x-shear can be expressed as a
fractional delay of one row at a time. In turn, a fractional delay of a ID signal can be
accomplished using a shift-invariant filter. Similar decompositions are known for 3D images
and shear operations.
Resampling a digital image can also be usually performed using lower dimensional
operations, which can often be shift invariant. For example, image upsampling along one
coordinate by an integer factor U may be decomposed into U different, computationally
efficient, fractional delays. This is essentially, the so-called polyphase decomposition, well-
known in digital signal processing. Rational, but non-integer resampling along one coordi-

nate can be decomposed into a cascade of integer down and upsarnpling, each of which is
efficiently performed.
More general digital image resampling can also be decomposed into lower dimen-
sional operations. Consider the resampling of a digital ?,D image / from a sampling pattern
with samples Si lying on a family of curves, denoted CFi, to another another pattern with
points S-2 lying on a different family of curves, CF2, producing the digital image h. If the two
families of curves intersect at sufficient density the method mat was described with reference
to Fig. 23 for resampling from one fan to another rotated fan can be used for general curves.
Otherwise a third family of curves CF3 can be introduced, which intersects both CFj and
CF2 at a desired density. The digital image can then be resampled from CFi to its intersec-
tions with CF3, then to the intersections of CF3 with CF^, and finaly on CF2 to the desired
sampling pattern.
This process generalizes to 3D, for example by considering surfaces instead of curves,
to first reduce the process to one of resampling on surfaces, and then using the 2D process to
further reduce resampling on surfaces to resampling on curves.
Digital coordinate transformations and resampling can often be combined, improving
the computational efficiency. For example, the digital x-shear-scale operations used in the
shear-scale algorithm shown in Fig. 13 can be decomposed into resampling operations on
individual horizontal lines.
1 Divergent-beam Fast Hierarchical Backprojection Algo-
rithms
1.1 Fan-beam projection and backprojection
Consider the case of equiangular-spaced detectors located on a circle around the ob-
ject. The fan-beam tomographic projection, at a source-angle /3, of a two-dimensional image
/(i, y) is denoted by (7^/) (7) and is defined as the set of line integrals along the rays of
a fan. parameterized by 7. centered at the source position on a circle of radius D from the
origin. The function f(x) is assumed to be zero-valued outside a disc of radius B.


The fan-beam projection at source-angle {} and fan-angle 7 is {TZpf)^) — f^Q f{V(/3, j,T))d,T.
Since / is zero outside the disc of radius R, the integral needs only be performed within the
disc i.e. between T$T and TEND-
In computed tomography with the fan-beam geometry, projections are available at a
set of discrete source angles {/3P : p — 1,2,..., P}, and within each fan the angles of the
rays are indexed by {jj : j — I, 2,..., J}. In the case of equiangular fan-beam geometry
[he detectors are equally distributed on the arc of a circle centered at the source, so the
fan-angles are evenly spaced. Ln the case of equispaced detectors, the detectors are equally
distributed on a line perpendicular to the line from the source position to the origin. The
fast backprojection algorithms for fan-beam geometry described here assume equiangular
distribution with J odd and 7^ = A7 ■ (j — (J -f l)/2). However, the algorithms may be
easily extended-to other fan-angle geometries.
The reconstruction algorithm from a set of P fan-beam projections {Jt^.(7) : i —
1, 2,..., Pj, may be expressed as the scaling and filtering of each fan-beam projection fol-
lowed by a weighted backprojection (5):

where W(T) is an appropriate weight function, T((x),/3)) and j((x),/3) are the distance
along the ray between source and image point £ for source position (3, and the fan angle of
that ray, respectively, and qp(j) are the weighted and filtered projections

where 5(7) is an appropriate filter,
Definition 1.1 The weighted backprojection of a function q at a single source angle (3 is
defined by

and the weighted backprojection operator is defined by

Therefore
It is easily shown that
where fC((3) denotes the rotation by /3 operator, and consequently that

Thus, as in the parallel-beam case, the weighted backprojection of P fan-beam projec-
tions may be expressed as the sum of weighted-zero-backprojected images as seen in Figure
4, with B0 — )%. In fact, close analogs of the backprojection operator defined in Equa-
tions (43) and (44) apply more generally to the reconstruction of functions in two and higher
dimensions from projections of a general form, with appropriate definition of the functions
7(z, f3) and T(x, p). The methods of the present invention extend to these other applications,
with appropriate definition of a coordinate transformation K.
1.2 Fast hierarchical backprojection
The Fast Hierarchical Backprojection Algorithm for the fan-beam geometry is similar
to that for the parallel-beam case. It combines the fan-beam projections in a ternary hierar-
chical structure, exploiting the fact that intermediate images formed by projections that are
close to each other in projection angle ft can be sampled sparsely. The block-diagram is
shown in Figure IS. Similar algorithms radix structures.
The equations that govern the combination of the underlying continuous images in

Figure 18 are as follows :

In the embodiment described next, it is assumed thai the source angles are uniformly
spaced in [0, 2ir) as follows :

The intermediate rotation angles are then chosen as in the parallel-beam case using
Equation (22), with A^ replacing Ag.
The constituent fans of an intermediate image in three levels of the algorithm with
block diagram Fig. 18 are illustrated in Fig. 20.
1.2.1 The fan-beam sampling scheme
One embodiment of the fan-beam algorithm us?s a sampling scheme of the intermedi-
. ate images in the algorithm derived by analogy to the parallel-beam case, which is described
next. Alternatives and improvements are described later.
A single-fanbeam-backprojection image (an image Wo? produced by weighted back-
projection of a single fanbeam ) has a structure amenable to sparse sampling. It follows from
Equation (43)that {Wpq)(V{p} 7, T)) = W(T)q(-y). A single sample at a particular T = 7',
on the ray indexed by 7 in a fan oriented at /? = 0 is merefore sufficient to specify the value
of the image along the entire ray: it is simply proportional to W{T) where 71 is the distance
from the source to the point in the ray.
An intermediate image in the algorithm is a sum of several rotatated versions of such
smgle-fanbeam-backprojection images, each generated from a constituent projection. We re-
fer to the fan at the angle that is at the center of the angular interval spanned by the constituent
projections as the central constituent fan. This may correspond to an actual projection - e.g..
in the case of a ternary algorithm, - or a to virtual projection - e.g., in the case of a binary

algorithm. For example, the fan in Fig. 20(a) is the central constituent fan for both Fig. 20(b)
and Fig. 20(c).
Recall that in the parallel-beam algorithm the intermediate images are sampled on
cartesian sampling patterns aligned with the central constituent. In the rotation-based parallel-
beam case the intermediate images are sampled along parallel vertical rays. The spacing of
samples along each of these parallel rays is chosen to guarantee that the other constituent
projections of the image are sufficiently sampled (by the Nyquist sampling criterion) along
that ray. The necessary spacing is exactly equal to the'spacing of the intersection points of
the vertical rays with a set of rotated parallel rays corresponding to the extremal constituent
projection, that is, the one farthest away in view angle from the center projection.
In direct analogy to the parallel-beam case, here the intermediate images are sampled
along the rays of the central constituent fan. The sampling points along the rays of the
central fan are the intersection points of the central fan with the extremal constituent fans.
This results in samples that are not uniformly spaced along each ray. Consider two fans
V(/3, jj, -)\Jj-i and V{/3', y'-, -)|/=i. Assuming the equiangular fanbeam geometry with an
odd number ,/ of detectors, y1 — A7 • (j — (J + l)/2). The points on the rth ray of the 8
-fan (corresponding to fan-angle yr) that intersect the /J'-fan, are determined by the equation
V{{3, yr, T) = V(P', y'j, T')|/=3, which has the solution

The function Tp^nT(~l') describes how the fan V(j3', -, -) varies along the rth ray of
the fan V{8, -, •)• It carries information about the varying sampling rate along a ray of the
fan. The steeper the slope, the sparser are the samples. The local sampling rate at any y1 is
proportional to (dT/dy')~l. A typical T(y') is displayed by the dashed line in Fig. 22.
Two modifications are made to Tp^>ilr (-/) to get a new sampling function T^p- lr (7'):
1. Intermediate images are to be sampled within the disc of radius R, with margins of
at least one sample on each ray outside the disc (one sample with T with T > TEND). Define y'END — TjfJ,i>rfr(T$ND) - the value of 7' corresponding to TBN>D.
Fans with adjacent source angles may lack intersection points with T > TEND- in order to
rectify that, the slope of T$fi>^T (7') for T > TEND is clamped to that at Tp$i^r {y'END)-

2. The final image is sampled on a Cartesian pattern with unit sampling intervals.
This is the smallest interval at which intermediate images need to be sampled. The point:
7^ at which a unit sampling race along the 7^ ray is achieved is determined by solving the
equation dI'(dj)Tp$*iT> (7) ~ 1/A7 for 7. To maintain the unit sampling interval constraint,
the slope of f'ptp'ilr (7') for 7' > 7^ is set to 1/A7.
Fig. 22 displays Tp^i^r(j) and Tptp>ilT(j) for a typical f3 and 7. It also displays the
T-values of the actual sampling points on the ray which are the set {Tpj*^ C?A7) : j c- Z}.
The integers j are chosen so that there are margins of at least one sample on each side of the
image boundary.
Clearly, the locations of the sample points prescribed by the principles outlined herein
can be computed from the geometry of the fans, Equation (21) for the set Ni>rr> of the con-
stituent projections, the selected rotation angles SijTn for the intermediate images (given in
Equation (22) for the case of equispaced view angles), and the chosen form for Tp^inr(j).
Separable rotation and up-interpolation for fan-beam sampling
As described in the Overview, the upsampling operations, and upsampling combined
with rotation operations can be decomposed into computationally efficient one dimensional
resampling operations.
1.2.2 Sampling Scheme Based On Local Fourier Structure
Given an image f(x) : Rn —> R, we want to find a sampling function t(x) : Wl —> R"
such that f(t(x)) has a small essential bandwidth and therefore can be sampled on the set of
points (£(m) : in shortly, knowing how f(x) is composed of its constituent projections, we find the matrix
function v(x) — Vvl][s) 111(f) t^at describes how the function / should be sampled locally
at the point, x. We then integrate this matrix function over the image domain to get the
sampling function t{x).
In our algorithms, we know that the intermediate image / is of the form : / =
;r;p/>C(Jp)y\Vip for some angles 5P assuming that the projections qp are sampled exactly at the Nyquisr rate, we know the lo-

cal structure of the spectral support at each point in the image /. Because of the fan-
backprojection of each band-limited projection qv, an image that is fan-backprojected from
a single projection has a negligible spatial bandwidth in the direction of the spoke of the fan,
while it's bandwidth in the perpendicular direction is inversely proportionate to the distance
from the vertex of the fan as seen in Figure 24. When a set of these fans are rotated and
added together, the local spectral support of this resulting image is the union of the spectral
supports of the individual fans.
For example when the constituent fans have source angles between /3min and /3mai, as
shown in Figure 25, the local spectral support at a point x in the image domain is a warped
wedge with radii oriented between 6min and B^x, and a radial bandwidth of QR/TQ at angle
6. Here Q, is the spatial bandwidth of the back-projected projection at the center of the unage
plane assuming the projections are sampled at the Nyquist rate. In the equiangular case Q is
number of samples of the projection divided by the length of the arc through the origin of the
image plane that the backprojected projection covers. Mathematically the spectral support is
where
Given this knowledge of the local two-dimensional fourier structure at a point x in
the image, the matrix function v(x) at that point is the sampling matrix that, if the spectral
support at that point were uniform across the whole image (such as in the parallel-beam
case), would efficiently sample the image. In the intermediate images of interest here, this
produces two distinct small-bandwidth and large-bandwidth sampling directions. We fix the
first coordinate to be the small-bandwidth direction.


Here Qmin £ Q#min and nmQ3. = Clo^. Angle 0mid refers to the angle of the projection
corresponding to the. source-angle (/3min + Pmax)/2, and Qmid. — &emid- These parameters,
chosen by geometrical arguments ©n the spectral support, ensure that the spectral support
upon transformation by the sampling matrix is restricted to [—7r, TT] X [—TT, n}.
The sampling function for the whole image is found by integrating this sampling ma-
trix function across the whole image — solving the set of differential equations ' £f- —
Vij(t.(n)) for i,j — 1,2. This may be solved numerically. Even if this exact prescribed
pattern is not used, the local Fourier support analysis will evaluate the effectiveness of any
sampling pattern.
The resulting sampling patterns for a few intermediate images are shown in Figs.
26(A) and 26(B). For the fan-beam case, this local Fourier-based method produces sampling
patterns that are similar to the ones resulting from the intersection-based methods described
earlier.d The previously described separable method to resample from one sampling pattern
to another can be used with these patterns also. This local Fourier sampling method may
be applied directly to find sampling schemes for arbitrary projection geometries over lines,
curves or planes over arbitrary dimensions. In the rase of the parallel-beam geometry, it
reduces to that discussed earlier.
1.2.3 Alternative Sampling Schemes
The sampling schemes in which the samples are located on a fan whose vertex is on the
source trajectory (called a "sampling fan") has some shortcomings. The sampling points are
chosen separately for each ray resulting in a scheme that is not necessarily optimal in a two
dimensional sense. Though the final image in the algorithm is sampled on a uniform rect-
angular grid with unit-spacing, the intermediate images using the above mentioned scheme
are sampled more densely than necessary in certain regions of the image (eg. close to the

vertex of the fan). In order to rectify this the above described scheme might be modified to
make sparse the sampling in such oversampled regions. Two such possibilities are illustrated
in Figure 27. Both these possibilies incorporate the fact while in the first level the image is
sampled efficiently on a fan, in the final level it needs to be sampled on a rectangular grid.
These schemes attempt to incorporate the gradual transition from sampling on a fan to a grid.
In Fig. 27 (A) samples are selected on a pseudo-fan whose vertex is located further
from the origin than the source radius. In successive levels intermediate images are formed
from fans from a larger range of source angles and the distance of the vertex of the pseudo-
fan from the origin increases. In the final l&vel the vertex is at infinity; i.e. the rays are the
parallel lines of the rectangular grid. In Fig. 27 (B) the samples are located, not on a fan,
but on a beam that is less divergent nearer the source position. In successive levels the beam
becomes more parallel and less divergent. These or other sampling schemes that take a two
dimensional (or, for-3D images, a 3D) point of view will contribute to a faster algorithm.
1.2.4 Other Optimizations
• la the fan-beam case with a full trajectory, source angles are distributed in fO, 2ir).
This allows us to take advantage of rotations by — TT/2, — n and —37r/2 that involve no
interpolation but only a rearranging of pixels.
• The positions of points in the sampling patterns for all levels of the algorithm can
be pre-computed and stored in a look up table, or computed on-the-fly along with the
processing of data. In either case, computation and/or storage can be reduced by taking
advantage of the smooth variation of sample position as a function of ray index, which
may be observed in Figs. 23(A) and 23(A).
1.2.5 Oversampled Fast Hierarchical Backprojection
As in the parallel-beam case, oversampling by a factor 7 accuracy. One way to achieve this in the fan-beam case, is to determine the denser sampling
patterns using a fan-angle spacing A'p = Ap - 7.

We Claim
1. A method (Fig. 5A) for creating a pixel image, f, from projections
( q1,..., qp ) comprising the steps of:
(a) producing (100) intermediate-images (I1,1..,I1,P) from selected
projections (q1,...,qp) ;
(b) performing digital image coordinate transformations (102) on selected
intermediate-images [ l11... ,I1,p), the parameters of coordinate transformations being
chosen to account for view-angles of the projections from which the intermediate images
have been produced, and for the Fourier characteristics of the intermediate-images;
(c) aggregating subsets of the transformed intermediate-images (104)
produced in step (b) to produce aggregate intermediate-images ( I2,1...,I2,p/2} '■> and
(d) repeating steps (b), and (c) in a recursive manner until all of the
projections and intermediate images have been processed and aggregated to form the
pixel image, f;
wherein the coordinate transformation parameters are chosen so that the
aggregates of the intermediate-images (104) may be represented with desirable accuracy
by sparse samples.
2. A method for creating a pixel image, f, from projections (ql,...,qp)
comprising the steps of:
(a) producing (99) a plurality of intermediate-images (I1,1,..., I1,p ), with at
least one corresponding to a non-Cartesian and/or non-periodic sampling pattern;
(b) performing digital image upsampling or downsampling (106) on selected
intermediate-images (I1,1...,11,p );

(c) performing digital image coordinate transformations on
upsampled/downsarnpled intermediate-images;
(d) aggregating (110) subsets of the transformed intermediate-images
produced in step (c) to produce aggregate intermediate-images (I2- . •,I2,P/2] ; and
(e) repeating steps (b), (c) and (d) in a recursive manner until all of the
projections and intermediate images have been processed and aggregated to form the
pixel image, f;
wherein at least one of the digital image coordinate transformations is performed
with a non-Cartesian and/or non-periodic sampling pattern, and the coordinate
transformation parameters are chosen so that the aggregates of the intermediate-images
may be represented with desirable accuracy by sparse samples.
3. The method as claimed in claim 1 or 2, wherein said aggregation is
performed by adding digital images.
4. The method as claimed in one of the preceding claims, wherein at least
one intermediate image is produced in step (a) by at least one of backprojection and
weighted backprojection (180,182) of selected projections.
5. The method as claimed in one of the preceding claims, wherein at least
one intermediate image is produced by at least one of backprojection and weighted
backprojection (180, 182) of two or more selected projections in step (a).
6. The method as claimed in one of the preceding claims, wherein at least
one aggregate intermediate image is formed by aggregating three or more selected
transformed intermediate images.
7. The method as claimed in one of the preceding claims, wherein the digital
image coordinate transformations are performed using digital filtering.

8. The method as claimed in one of the preceding claims, wherein selected
coordinate transformations are digital image rotations.
9. The method as claimed in one of the preceding claims, wherein selected
coordinate transformations are digital image shearing (120, 122), or shear-scaling.
10. The method as claimed in one of the preceding claims, wherein selected
coordinate transformations are at least one of upsampling (101, 106) and downsampling
(109) of the digital images.
11. The method as claimed in one of claims 7 to 10, wherein at least some of
said digital filtering is performed by one-dimensional digital filters.

12. The method as claimed in one of claims 7 to 11 in which at least some of
said digital filtering is performed by shift-invariant digital filters.
13. The method as claimed in one of claims 7 to 12, wherein at least some of
said digital filtering is recursive.
14. The method as claimed in one of claims 7 to 13, wherein at least some of
said digital filtering is implemented using a fast Fourier transform (FFT).
15. The method as claimed in one of the preceding claims, wherein a selected
degree of oversampling is applied to at least one of selected intermediate images,
transformed intermediate images, and aggregate intermediate images;
16. The method as claimed in one of the preceding claims, in which non-
Cartesian sampling patterns are used.

17. The method as claimed in one of the preceding claims, wherein at least
one of formation of selected intermediate images in step (a), selected coordinate
transformations, and aggregation steps are combined within a level, or across adjacent
levels of the hierarchy or recursion.
18. The method as claimed in one of the preceding claims, wherein at least
one intermediate image is at least one of weighted before and after performing digital
image coordinate transformations.
19. A method for creating a pixel image f from projections (ql,...,qp)
along a collection of lines, curves, or surfaces comprising the steps of:

(a) producing (184) intermediate images (I1,1..., Il,P);
(b) performing digital image resampling on selected intermediate images
(186), the location of samples being chosen to account for the view-angles of the selected
projections and for the Fourier characteristics of the intermediate images,
(c) aggregating (190) selected subsets of the resampled intermediate-images
to produce aggregate intermediate-images [ I2,1, •......I2,p/2) ; and
(d) repeating steps (b) and (c) in a recursive manner, at each level of the
recursion increasing the density of samples of the intermediate images, until all of the
projections and intermediate images have been processed and aggregated to form the
pixel image;
wherein the sampling scheme is chosen so that aggregates of the resampled
intermediate-images may be represented with desirable accuracy by sparse samples.
20. The method as claimed in claim 19, wherein at least one intermediate
image is produced in step (a) by weighted backprojection (180,182) of selected
projections.

21. The method as claimed in claim 19, wherein at least one intermediate
image is formed by weighted backprojection (180,182) of two or more selected
projections in step (a).
22. The method as claimed in claim 19, wherein at least one aggregate
intermediate image is formed by aggregating three or more selected transformed
intermediate images in step (c) or (d).
23. The method as claimed in one of claims 19 to 22, wherein the intermediate
images have samples that lie on a family of lines, curves, or surfaces.
24. The method as claimed in one of claims 19 to 23, wherein the digital
image resampling is performed by a sequence of lower-dimensional digital filtering
operations by utilizing intermediate sampling schemes that lie on the intersections of the
families of lines, curves or planes.

25. The method as claimed in one of claims 19 to 24, wherein a selected
degree of oversampling is applied to at least one of the selected resampled intermediate
images and aggregated intermediate images.
26. The method as claimed in one of claims 19 to 25, wherein said
aggregation is performed by adding digital images.
27. The method as claimed in one of claims 19 to 26, wherein at least one of
formation of selected intermediate images in step (a), resampling, and aggregation steps
are combined within a level, or across adjacent levels in the hierarchy or recursion.
28. The method as claimed in one of claims 19 to 27, wherein at least one
intermediate image is at least one of weighted before and after resampling step (b).

29. The method as claimed in one of claims 19 to 28, wherein changes in
sampling density are accomplished by digital filtering.


ABSTRACT

METHODS FOR CREATING A PIXEL IMAGE FROM PROJECTIONS
In the present invention pixel images (116) are created from projections by
backprojecting selected projections to produce intermediate images, and performing
digital image coordinate transformations (102) and/or resampling on selected
intermediate images. The digital image coordinate transformations (102) are chosen to
account for view angles of the constituent projections of the intermediate images and for
their Fourier characteristics, so that the intermediate images may be accurately
represented by sparse samples. The resulting intermediate images are aggregated into
subsets (104), and this process is repeated in a recursive manner until sufficient
projections and intermediate images have been processed and aggregated to form the
pixel image (116). Digital image coordinate transformation can include rotation (102),
shearing, stretching, contractions, and the like. Resampling can include up-sampling,
down-sampling, and the like.

Documents:

00725-kolnp-2006-abstract.pdf

00725-kolnp-2006-claims.pdf

00725-kolnp-2006-description complete.pdf

00725-kolnp-2006-drawings.pdf

00725-kolnp-2006-form 1.pdf

00725-kolnp-2006-form 3.pdf

00725-kolnp-2006-form 5.pdf

00725-kolnp-2006-international publication.pdf

00725-kolnp-2006-international search report.pdf

00725-kolnp-2006-pct request form.pdf

00725-kolnp-2006-priority document.pdf

725-KOLNP-2006-(07-03-2012)CORRESPONDENCE.pdf

725-KOLNP-2006-ABSTRACT.pdf

725-KOLNP-2006-AMANDED CLAIMS.pdf

725-KOLNP-2006-ANNEXURE FORM 3.pdf

725-kolnp-2006-assignment.pdf

725-KOLNP-2006-CORRESPONDENCE 1.1.pdf

725-KOLNP-2006-CORRESPONDENCE 1.2.pdf

725-KOLNP-2006-CORRESPONDENCE 1.4.pdf

725-KOLNP-2006-CORRESPONDENCE-1.3.pdf

725-kolnp-2006-correspondence-1.5.pdf

725-KOLNP-2006-CORRESPONDENCE.pdf

725-KOLNP-2006-DESCRIPTION (COMPLETE).pdf

725-KOLNP-2006-DRAWINGS.pdf

725-KOLNP-2006-EXAMINATION REPORT.pdf

725-KOLNP-2006-FORM 1.pdf

725-kolnp-2006-form 18-1.1.pdf

725-KOLNP-2006-FORM 2.pdf

725-KOLNP-2006-FORM 5.pdf

725-kolnp-2006-granted-abstract.pdf

725-kolnp-2006-granted-claims.pdf

725-KOLNP-2006-GRANTED-DESCRIPTION (COMPLETE).pdf

725-kolnp-2006-granted-drawings.pdf

725-kolnp-2006-granted-form 1.pdf

725-kolnp-2006-granted-form 2.pdf

725-KOLNP-2006-GRANTED-SPECIFICATION.pdf

725-kolnp-2006-others-1.1.pdf

725-KOLNP-2006-OTHERS.pdf

725-KOLNP-2006-PA.pdf

725-KOLNP-2006-PETITION UNDER RULE 137.pdf

725-KOLNP-2006-REPLY TO EXAMINATION REPORT.pdf


Patent Number 252715
Indian Patent Application Number 725/KOLNP/2006
PG Journal Number 22/2012
Publication Date 01-Jun-2012
Grant Date 28-May-2012
Date of Filing 27-Mar-2006
Name of Patentee THE BOARD OF THE TRUSTEES OF THE UNIVERSITY OF ILLINOIS
Applicant Address 506 SOUTH WRIGHT STREET, 352 HENRY ADMINISTRATION BUILDING, URBANA IL
Inventors:
# Inventor's Name Inventor's Address
1 GEORGE, ASHVIN, K 1001 W. CLARK STREET APPT C1, URBANA, IL 61801
2 BRESLER. YORAM 414 BROOKENS DRIVE, URBANA, IL 61801
PCT International Classification Number G06T
PCT International Application Number PCT/US2004/029857
PCT International Filing date 2004-09-09
PCT Conventions:
# PCT Application Number Date of Convention Priority Country
1 60/501,350 2003-09-09 U.S.A.