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 upsampling, downsampling, and the like. 
Full Text  This is a continuationinpart 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 wellknown technique underlying nearly all of the di agnostic imaging modalities including xray 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) parallelbeam, in which the lineintegrals are performed along sets of parallel lines; (ii) divergentbeam, in which the lineIntegrals 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 lineintegral projections. Another problem of tomographic reconstruction is to reconstruct a 3D image from a set of its surfaceintegral 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 parallelbeam 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 Npixel image in 2D, and at least as NA for an N x. N x Nvaxel image in 3D. Thus, doubling the image resolution from N to 2N results in roughly an 8fold (or 16fold, 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 multirow 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 conebeam 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 parallelbeam tomography using a multilevel decomposition, producing at each stage an image covering the entire fieldofview, 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 upsampling, downsampling and the like. Projections can be created from a pixel image by performing digital image coordinate transformation and/or resampling and/or decimation and reprojecting 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. 20A20C; Figs. 23A, 23B, 23C and 23D illustrate sampling points for rotation for upsampling 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 twovariable functions are denoted variously, depending on conve nience. The following three notations of function f are equivalent: f(x1,x2), f (f), and f(£]). Continuousdomain and discretedomain 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 2D 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 lineintegral projections, e.g., with a divergent beam geometry, 14 to a projection preprocessor 16. The projection preprocessor 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 preproccssor 16 is a collection of preprocessed 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 beamhardening, or as part of a multistep 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 multistep 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 multidimensional signal processing, for example as described in D. Dudgeon and R, Mersereau, Multidimen sional Digital Signal Processing. Englewood Cliffs: PrenticeHall, 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. 2A2C and 3A3C, 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 zeroth 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. Upsampling and downsampling are special cases of resampling on a denser or sparser sampling pattern, respectively. Further, upsampling or downsampling by a factor of 1 involves no change in the digital image, is considered a form of upsampling or downsampling. Digital image addition refers to to a pointbypoint 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, Convolutionbased interpolation for fast, highquality rotation of images. IEEE Transactions Image Processing, Vol. 4, pp. 13711381, 1995. Some digital image coordinate transformations are illustrated in Figs. 2A2C 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. 2B and 2C 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. 2A. 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 upsampling or downsampling. 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 discretetime 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 unweighted 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 rotationbased 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,1I1,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 wedgeshaped 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 Fourierdomain 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 brokenline rectangle) of the intermediate images can be minimized, allowing sparse sampling and reducing compu tational requirements. Fig. 7(C) shows the spacedomain 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 discretetime Fourier transforms (DTFTs) of the digital images Ilm. Fig. 10(A) describes an algorithm for another embodiment of the present invention, using ternary twoshearbased 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 shearscale backprojection, and Fig. 1 IB illustrates an equivalent hierarchical shearscale backprojection. In Fig. 11A, the plurality of projections q1:..., q4, are backprojected at 140 and subjected to a shearscale 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 shearscale 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 shearscale backprojection. Projections qe1, ...,q0l are backprojected at 160 and subjected to a shearscale digital image coordinate transformation at 162. Subsets of the resulting images are aggregated at 164, and the resulting intermediate images are sub jected to upsampling 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 twoshear 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, xcoordinate shear comprising the twoshear digital image coordinate transformation 108 (shown in Fig. 10B) producing a shearscale 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 shearscale 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 xshear step 168 Sx',m, the two can be combined for computational effi ciency into a single digital image shearscale. Noncartesian 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 noncartesian grid, Noncartesian 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 nonCartesian sampling patterns can be executed ef ficiently using one dimensional shiftinvariant 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 nonCartesian 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 Xray transform that arises for example in Positron Emission Tomography, and the 3D Radon transform, which arises in magnetic resonance imaging. The threedimensional (3D) Xray transform is a collection of integrals along sets of parallel lines, at various orientations, in 3D. Each 3D Xray projection is a twodimensional function that can be characterized by the 3D angle at which the lines are oriented. The block diagram for hierarchical backprojection for 3D Xray data is similar to the ones previously described, such as Figures 8. 10, 13, 15 or 16. The intermediate images in this case are threedimensional, not twodimensional 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 sparsesampling (slow) direction increases to accommodate the increasing bandwidth in that direction as explained by the Fourier analysis of 3D Xray 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 onedimensional operations, such as shears and shearscales, 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 onedimensional function — a collection of in tegrals along sets of parallel 2D planes, parameterized by the displacement of the plane from the origin. The viewangle of the projection is that of the 3D orientation angle of a vector perpendicular to the set of planes. The blockdiagram 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 viewangles 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 Bspline transforms for continuous image represen tation and interpolation, IEEE Transactions on Pattern Analysis and Machine Intelligence, Vol. 13, pp. 277285, 1991; M. Unser, A. Aldroubi, and M. Eden, Bspline signal, process ing: Parr Ithf.ory, IEEE Transactions Signal Processing, Vol. 41, pp. 821832, 1993; and M Unser, A. Aldroubi, and M. Eden, Bspline signal processing: Part II  efficient design and application'!, IEEE Transactions Signal Processing, Vol. 41, pp. 834848, 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 fanbeam and conebeam, with arbitrary source trajectories. The common fanbeam 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 fanbeam 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 imagedisc 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 rotationbased 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 onebeforeiast 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. 20A20C 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 intersectionbased 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 backprojection, 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 Fourieranalysis 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 backprojections on arbitrary trajectories. in both two and three dimensions. For the case of backprojection on nonperiodic 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 bowtie 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 Fourierbased method for determining the sampling patterns for resampling and coordinate transformation for hierarchical fanbeam 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 DivergentBeam Algorithms The methods described for fanbeam backprojecticn, extend directly to other divergent beam geometries, including 3D conebeam, which is one of the most important in modern diagnostic radiology as will be described. Similarly to the fanbeam geometry, in which a source of divergent rays moves on a trajectory (i.e., a circle) around the imaged object, in the general divergentbeam 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 divergentbeam projections are zerobackprojected at 99 with the appropriate divergentbeam singleview 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 divergentbeams of the intermediate images may not simply rotated as in the fanbeam 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 rearrangement which replaces a n/2 rotation in the fanbeam algorithm. As in the fanbeam case, the. initial intermediate images are processed hierarchically by the algorithm. Analogously to the fanbeam 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 divergentbeams 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 fanbeam 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 fanbeam case. This matrix function is then integrated (possibly numerically) over the image domain to determine the position of the samples. A separable method of upsampling combined with rotation and translation onto a new sampling beam is achieved similarly to the fanbeam case. It reduces the 3D coordinate transformation and resampling operations to a sequence of ID resampling operations. As shown in Figure 29 (a), each divergentbeam may be regarded as the intersection of a set o( vertical planar fanbeams that are distributed in azimuthal angle, with a set of tilted planar fanbeams at different elevation angles . The steps of the separable coordinate transformation are as follows (as shown in Figure 29 (b)) 1. The original divergentbeam is resampled onto the set of intersection of the rays of the original with the vertical planes of the new divergentbeam . 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 fanb&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 divergentbeam can be expected to achieve large speedups with desirable accuracy. As in the fanbeam case one might use a pseudobeam 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 resamplingbased 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 rotationbased backprojection 4, but differs from it in two respects. First, in the first step the pth (possibly processed) projection is subjected to a weightedbackprojection 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 sourceangle. 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 singleprojection intermediate images is chosen to be an efficient and sparse pattern, so will usually be different for each projec tion, and often nonCartesian. 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 resamplingbased 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 nonhierarchical 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 nonCartesian 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 fanbeam 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 nonperiodic 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 blockdiagram 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 unrotated 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 (lowpass 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 viewangles, 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 onedimensional, signal. These onedimensional 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 blockdiagram is a flow transposition of a backprojection blockdiagram, every branch of the reprojection blockdiagram has a corre sponding branch in the backprojection blockdiagram, 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 twoshear backprojection algorithm in Fig 10, the coordinate transformations in the second level (202) is an xshear followed by a yshear, and the coordinate transformation in the last level (206) is an xshear,fo!iowed by a yshear, and a fractional decimation in x (These three operations can be reduced to a shearscale). 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 shearseal^ version of this algorithm, (corresponding to the shearscale 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 shearscales. The shearscale used in the reprojection algorithm is the mathematical adjoint of the corresponding shearscale 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 blockdiagram of a decimationbased weighted reprojection. It it the flowgraph 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 onedimensional decimations (lowpass 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 onedimensional 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 SheppLogan 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 rootmean  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 parallelbeam data, the test image was of size 256 x 256, the number of view an gles was 486, and the SheppLogan 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 twoshear algorithm is represented by the circles and the shearscale algorithm by the squares. Each algorithm is run using two types of filters — a thirdorder (dashed line) and fifthorder (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 nonoversampled 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 parallelbeam data, are displayed in Fig. 35. Columnwise from left to right, they are output images from the con ventional backprojection, the twoshear and the shearscale 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 fanbeam geometry was successfully tested on the 512 x 512 2D shepplogan 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 fanbeam algorithm that used linear interpolation. In all the experiments the projections were, filtered with the SheppLoga.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 xaxis of the psfs of the conventional algorithm, the nonoversampled 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 intersectionbased 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 (datadependant) 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 unenhanced 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 parallelbeam 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 socalled ramp or modified ramp filter and then backprojecting those filtered projections. Additional preprocessing 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 lengthP vector 0), is described by: where Definition 0.1 the backprojection operator is defined by The 0(N2logP) coordinatetransformationbased 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) rotationbasec hierarchical backprojection algorithms of the present invention are based on decomposition of backprojection in terms of the rotation of images. The rotationby0 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 viewangle 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 viewangles, 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 onedimensional Fourier transform (J^i) of a projection Vgf with the twodimensional Fourier transform (T2) of the image /. The projectionslice theo rem 4 says that The backprojection algorithm of Fig, 4 can be interpreted in the Fourier domain. The backprojectionatzero 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 onesided 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 viewangles 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 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 rotationbased 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 xi coordmate will be needed for initial intermediate images formed, as in Fig 3b, by backprojcciing more than one projection. The sum of two virtual intermediate images (lijin = Iii2m\ + hx^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 Fourierdomain support of / and /. TFig. 7(C) shows the spacedomain 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(N2'/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 rotationbased 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 tinrotated one. This partic ular combination of binary and ternary stages results in a set of viewangles 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 viewangles 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 zeroangle 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 upsampling 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 viewangles 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 slowsampling 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 backprojectedatzero (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 £(,3mi = 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 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 xshear andyshear operators are defined by (Sx(a)f)(x) = f(S:x) iSv{a)f){x) = f(S«x) where the xshear andyshear matrices are S°  [\ f ] and S^  [£ §], respectively Tliis alternative embodiment is derived from the twoshear 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 twoshear digital image coordinate transforma tion, the twicesheared 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 noninteger resampling every time a twoshear 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 singleview 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 twoshear 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 twoshear 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. Shearscale algorithms The shearscalebased hierarchical backprojection algorithm is based on the definition of backprojection as the sum of singleprojection images that are scaled (stretched/contracted) and sheared. Definition 0,8 The xshearscale operator C(a) is defined by where the xshearscale 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 shearscales is also a shearscale. Using this property, Equation (23) can be rewritten in a hierarchical algorithmic structure. Fig. 11(A) is the blockdiagram 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 projectionangles 8P and the unknown shearscale parameters 5^, whose solution always exists. It is computationally beneficial (in terms of operations or storage space) to use an integer scalefactor rather than a noninteger one. For example, it is beneficial to use x h+l,m — 'SI(tti3m2)i'[3m2 + //,3m1 + 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 Similarly to the derivation of Equation (12), it can be shown that the intermediate shearfactors depend on the {PitTn} as follows : The freedom in the choice of the shear coefficients ai 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^mx 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 bandedges 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 spectralsupport .endpoints are denoted as Ef+l m,E}+1 , and Ef_l:m i.e.. Ell+1 m = Ei:imi. 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 shearfactors projection algorithm (Fig. 13). The inputs to the algorithm are the sets of endpoints 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 yupsampling 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 lowvalued denominators for the same reason. The problem of choosing these factors is slightly simpler in the case of the shear scale algorithm than the twoshear 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 Ylfit £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 ybandwidth 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 twoshear rotation algorithm, the problem is complicated slighdy by the fact that the use of twoshear rotation causes a defacto fractional up and down sampling of intermediate images. In particular, the slowsampling 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 ybandwidth 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 slowdirection 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 twoshear 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 twoshear algorithm therefore involves only one additional level of fractional resampling. This x coordinate resampling is combined with the digital 2shear in the Lth level for improved computational efficiency. The block diagram of this algorithm is in Fig. 15. In the shearscale hierarchical backprojection algorithm, the xupsampling operations are combined with the ^shearscales (C(a) at the beginning of the algorithm) to avoid an extra level of resampling. It is easy to see that :cupsampling is simply a special case of xshearscaling, with the shear factor set to 0. The combination of these two operations is still an xshearscale: C[a\, cr2)C(l/t/, 0) = C{a\/U, cr2/U). The downsampling at: the end of the algorithm is combined with the sshears in Lth level ( by a zdownsampling by U is effectively an xshearscale: 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 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 fastsampling 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 firstlevel images, B0qsp,p = 1,.... P. The block diagram of this ternary oversampled shearscalebased 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 xupsampling/downsampling, or xshearscales with xupsampling/downsampling. The nonoversampled twoshear algorithm of Fig. 10 involves a ysheai 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 xshear of this stage — into a single one. This leaves the length of the cascade of interpolations unchanged from the nonoversampling case. Hierarchies of Arbitrary Radixes/Branching factors All these algorithms can be easily modified for the case when the set of viewangles 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 viewangles, 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 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 righthand side of Equation (34). 0.1 Hierarchical Algorithms Based On Other Image Transformations The backprojection 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 inversetranspose 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 singleaxis 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 loworder 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 yshear of the digital image can be achieved by fractional delays of the individual vertical lines in a digital image array. Likewise, a digital xshear 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 shiftinvariant 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 socalled polyphase decomposition, well known in digital signal processing. Rational, but noninteger 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 S2 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 xshearscale operations used in the shearscale algorithm shown in Fig. 13 can be decomposed into resampling operations on individual horizontal lines. 1 Divergentbeam Fast Hierarchical Backprojection Algo rithms 1.1 Fanbeam projection and backprojection Consider the case of equiangularspaced detectors located on a circle around the ob ject. The fanbeam tomographic projection, at a sourceangle /3, of a twodimensional 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 zerovalued outside a disc of radius B. The fanbeam projection at sourceangle {} and fanangle 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 fanbeam 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 fanbeam geometry [he detectors are equally distributed on the arc of a circle centered at the source, so the fanangles 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 fanbeam geometry described here assume equiangular distribution with J odd and 7^ = A7 ■ (j — (J f l)/2). However, the algorithms may be easily extendedto other fanangle geometries. The reconstruction algorithm from a set of P fanbeam projections {Jt^.(7) : i — 1, 2,..., Pj, may be expressed as the scaling and filtering of each fanbeam 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 parallelbeam case, the weighted backprojection of P fanbeam projec tions may be expressed as the sum of weightedzerobackprojected 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 fanbeam geometry is similar to that for the parallelbeam case. It combines the fanbeam 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 blockdiagram is shown in Figure IS. Similar algorithms 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 parallelbeam 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 fanbeam sampling scheme One embodiment of the fanbeam algorithm us?s a sampling scheme of the intermedi . ate images in the algorithm derived by analogy to the parallelbeam case, which is described next. Alternatives and improvements are described later. A singlefanbeambackprojection 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 smglefanbeambackprojection 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 parallelbeam algorithm the intermediate images are sampled on cartesian sampling patterns aligned with the central constituent. In the rotationbased 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 parallelbeam 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, )\Jji 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 fanangle 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 Tvalues 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 upinterpolation for fanbeam 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 cal structure of the spectral support at each point in the image /. Because of the fan backprojection of each bandlimited projection qv, an image that is fanbackprojected 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 backprojected 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 twodimensional 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 parallelbeam case), would efficiently sample the image. In the intermediate images of interest here, this produces two distinct smallbandwidth and largebandwidth sampling directions. We fix the first coordinate to be the smallbandwidth direction. Here Qmin £ Q#min and nmQ3. = Clo^. Angle 0mid refers to the angle of the projection corresponding to the. sourceangle (/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 fanbeam case, this local Fourierbased method produces sampling patterns that are similar to the ones resulting from the intersectionbased 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 parallelbeam 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 unitspacing, 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 pseudofan 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, for3D images, a 3D) point of view will contribute to a faster algorithm. 1.2.4 Other Optimizations • la the fanbeam 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 precomputed and stored in a look up table, or computed onthefly 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 parallelbeam case, oversampling by a factor 7 accuracy. One way to achieve this in the fanbeam case, is to determine the denser sampling patterns using a fanangle 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) intermediateimages (I1,1..,I1,P) from selected projections (q1,...,qp) ; (b) performing digital image coordinate transformations (102) on selected intermediateimages [ l11... ,I1,p), the parameters of coordinate transformations being chosen to account for viewangles of the projections from which the intermediate images have been produced, and for the Fourier characteristics of the intermediateimages; (c) aggregating subsets of the transformed intermediateimages (104) produced in step (b) to produce aggregate intermediateimages ( 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 intermediateimages (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 intermediateimages (I1,1,..., I1,p ), with at least one corresponding to a nonCartesian and/or nonperiodic sampling pattern; (b) performing digital image upsampling or downsampling (106) on selected intermediateimages (I1,1...,11,p ); (c) performing digital image coordinate transformations on upsampled/downsarnpled intermediateimages; (d) aggregating (110) subsets of the transformed intermediateimages produced in step (c) to produce aggregate intermediateimages (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 nonCartesian and/or nonperiodic sampling pattern, and the coordinate transformation parameters are chosen so that the aggregates of the intermediateimages 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 shearscaling. 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 onedimensional 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 shiftinvariant 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 viewangles of the selected projections and for the Fourier characteristics of the intermediate images, (c) aggregating (190) selected subsets of the resampled intermediateimages to produce aggregate intermediateimages [ 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 intermediateimages 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 lowerdimensional 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 upsampling, downsampling, and the like. 

00725kolnp2006description complete.pdf
00725kolnp2006international publication.pdf
00725kolnp2006international search report.pdf
00725kolnp2006pct request form.pdf
00725kolnp2006priority document.pdf
725KOLNP2006(07032012)CORRESPONDENCE.pdf
725KOLNP2006AMANDED CLAIMS.pdf
725KOLNP2006ANNEXURE FORM 3.pdf
725KOLNP2006CORRESPONDENCE 1.1.pdf
725KOLNP2006CORRESPONDENCE 1.2.pdf
725KOLNP2006CORRESPONDENCE 1.4.pdf
725KOLNP2006CORRESPONDENCE1.3.pdf
725kolnp2006correspondence1.5.pdf
725KOLNP2006CORRESPONDENCE.pdf
725KOLNP2006DESCRIPTION (COMPLETE).pdf
725KOLNP2006EXAMINATION REPORT.pdf
725kolnp2006form 181.1.pdf
725kolnp2006grantedabstract.pdf
725kolnp2006grantedclaims.pdf
725KOLNP2006GRANTEDDESCRIPTION (COMPLETE).pdf
725kolnp2006granteddrawings.pdf
725kolnp2006grantedform 1.pdf
725kolnp2006grantedform 2.pdf
725KOLNP2006GRANTEDSPECIFICATION.pdf
725KOLNP2006PETITION UNDER RULE 137.pdf
725KOLNP2006REPLY TO EXAMINATION REPORT.pdf
Patent Number  252715  

Indian Patent Application Number  725/KOLNP/2006  
PG Journal Number  22/2012  
Publication Date  01Jun2012  
Grant Date  28May2012  
Date of Filing  27Mar2006  
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:


PCT International Classification Number  G06T  
PCT International Application Number  PCT/US2004/029857  
PCT International Filing date  20040909  
PCT Conventions:
