CA2324832A1

CA2324832A1 - Real-time image reconstruction for computed tomography systems

Info

Publication number: CA2324832A1
Application number: CA002324832A
Authority: CA
Inventors: Roger Lecomte; Germain Leger; Vitali Selivanov
Original assignee: Universite de Sherbrooke
Priority date: 2000-10-31
Filing date: 2000-10-31
Publication date: 2002-04-30

Abstract

The present invention discloses a method and a system for real-time reconstruction of tomographic images from projection data-generated by computed tomography scanners. The method is based on the theory of pseudo-inverse matrices applied to tomographic image reconstruction, and it allows immediate, projection-by-projection or event-by-event updating of the reconstructed tomographic image at the moment the measuring instrument supplies additional individual projection data. The data processing involves simple arithmetic, which can be implemented into general-purpose computers in the first case, and dedicated programmable logic or application specific circuits in the second case. Models of the underlying physical processes, such as the photon emission and detection processes or the spatially variant system response, can be easily included into image reconstruction through the system matrix. Corrections for other physical and instrumental factors affecting the accuracy of projection data, such as attenuation in the body or detector efficiency normalization, can be merged with the system matrix beforehand or be performed "on the fly" by simple multiplication of the matrix elements.
Real-time updating of the reconstructed tomographic image can be performed as soon as measured projection data become available. No storage of the measured raw data in the form of list mode events or sinogram is required, and the user of the tomographic system may choose to store reconstructed images only, thus achieving optimum compression of the tomographic data.

Description

REAL-TIME IMAGE RECONSTRUCTION FOR
COMPUTED TOMOGRAPHY SYSTEMS
BACKGROUND OF THE INVENTION
1. Field of the invention:
The present invention relates to a method for conducting real-time matrix pseudo-inverse tomographic reconstruction.

2. Brief description of the prior art:
Tomography refers to the cross-sectional imaging of an object from either transmission, emission or reflection data collected from many different directions.
Tomographic imaging deals with reconstructing an image from such data, commonly called projections. From a purely mathematical standpoint, a projection at a given angle is the integral of the image in the direction specified by that angle. The solution to the problem of reconstructing a function from its projections dates back to the paper by Radon in 1917 [1], thus the denomination inverse-Radon transform. However, it is only in the early 1970's, with the invention of the X-ray computed tomographic (CT) scanner [2,3] and the development of reconstruction algorithms [4-7] that tomographic imaging came into widespread practice. Tomographic methods now find applications in a vast number of fields, including radio astronomy, seismology, nondestructive analyses, electron microscopy, and above all medical imaging. In fact, most of the powerful new medical imaging modalities that have been introduced during the last three decades, such as X-ray CT, Single Photon Emission Computed Tomography (SPELT), Positron Emission Tomography (PET), Magnetic Resonance Imaging (MRI) and 3D ultrasound (US), were the result of the application of tomographic principles.
A. Review of current tomographic reconstruction algorithms Comprehensive description of the existing tomographic reconstruction algorithms can be found in textbooks by Herman [8], Dean [9], Natterer [10], Kak and Slaney [11] and Fessler and Ollinger [12], and review articles by Vardi et al.
[13] and Ollinger and Fessler [14]. There are two main classes of tomographic reconstruction methods currently in use today: analytical reconstruction algorithms, which aim at solving the inverse-Radon transform, and iterative algorithms, which attempt to reach a solution by successive estimation of the underlying image and improving image estimates at every iteration.
A.1 Filtered Backprojection The most common tomographic reconstruction method nowadays is the filtered backprojection (FBP) algorithm, which makes use of the projection-slice theorem [11]. The algorithm proceeds by first transferring each measured projection into Fourier space, filtering with an appropriate filter in frequency space, transferring the result back into the projection space and then backprojecting the result onto the image grid (see Figure 1 ):
z=Back-Project 1-I~c~f~xF l~b~)~ (1) where b~ is a j-th projection vector (j=1,...,A), A is a total number of projection angles, F 1 is a 1-D Fourier transform, F ~ -~ is a 1-D inverse Fourier transform, c(f~ is a 1-D filter function in the Fourier space, commonly referred to as a "ramp" filter, and x = {xk : k =1,..., M~ is the image that has to be found, M
being the total number of image pixels (or voxels).
The FBP algorithm produces an image that is a linear combination of the projection data. It can be efficiently implemented in general-purpose computers, its memory requirements are modest and the computation is rather fast. It is the reconstruction method of choice for virtually all X-ray CT and most emission tomography (SPELT, PET) systems currently in operation. Since the FBP
algorithm is linear, it can be reduced to event-by-event reconstruction based on the principle of superposition, by deriving a filter function representing the contribution of an individual event to the image. However, its widespread popularity stems more from historical reasons of computational simplicity than any widely accepted advantage in image quality. In fact, there are several problems with this algorithm:
- data must be assumed to be uniformly distributed on projections (which is generally not the case with the cylindrical geometry), or must be rebinned with equal spacing to accommodate the convolution operation (performed in Fourier space);
- models for the detector response must be space invariant and can only be incorporated into the algorithm as a deconvolution with the attendant noise amplification;
- the ideal ramp filter amplifies high-frequency noise in the projection, and it must often be apodized with a window function to reduce noise in the reconstructed image;
- even though the intensity is known to be non-negative, the algorithm yields negative values, particularly if the data are noisy;

- streak artefact is present; and - event-by-event reconstruction involves some elaborate computation to align, rotate and scale the single event filter before summation to the image matrix.
A.2 Algebraic Reconstruction There are other tomographic image reconstruction methods called algebraic reconstruction techniques [15,16]. They aim at solving the system of linear equations P11x1 +P12x2 +...+ plMxM = bl P21x1 +P22x2 +...+ p2MxM = b2 (2) PNlxl + PN2x2 +... ~- pNMxM = bN
where b; : i =1,..., N are the measured projections (sinogram), z = ~xk : k =1,..., M} is the image that has to be found, N and M are the total number of tubes-of-response and image pixels (or voxels), respectively.
These methods are usually iterative in nature and tend to converge to the minimum norm least-squares solution as the iteration number goes to infinity [16]. Some of the drawbacks of the iterative algebraic methods are:
- highly computer intensive, large memory requirements;
- because the algorithms are inherently iterative, fast reconstruction (e.g., parallel, distributed computation) methods can hardly be implemented;
- the need to optimize iteration number and the order in which projections are being utilized in the reconstruction;
- the stochastic nature of projection data in emission tomography is not taken into account; and - the complete measured projection data set is required before starting reconstruction.
A.3 Statistical Reconstruction Statistical image reconstruction methods were introduced in an attempt to overcome some of the drawbacks of the previous methods, in particular taking into account the stochastic nature of the measured data in emission tomography [13,17] and the physical modeling of the detection process [18]. In general, statistical iterative reconstruction methods yield images with much superior visual quality than FBP reconstructed images, especially in the case of low projection statistics (see Figure 2). The elimination of negative values and proper modeling of the detector response avoid several of the artifacts present in FBP
images. Unfortunately, this superior image quality does not necessarily translate into better quantitation accuracy, as this was recently demonstrated for the pharmacokinetic modeling of dynamic image sequences [19]. The use of statistical reconstruction methods also raises other practical problems that still hamper its utilization outside the research environment. In addition to most of the drawbacks that are common with the iterative algebraic methods, the statistical methods have the following problems:
- in the case of unconstrained iterative estimation, the absence of reliable, operator-independent, criteria to stop iterations: if iteration number is too low, the spatial resolution is sub-optimal; if iteration number is too high, the spatial resolution is artificially enhanced at the expense of unacceptable noise amplification; and - in the case of so called Bayesian approach, it is necessary to optimize an extra parameter, called prior, that controls noise [20].
SUMMARY OF THE INVENTION
In an attempt to overcome the above discussed drawbacks of the prior art, there is provided, according to the present invention, a method for conducting real-time matrix pseudo-inverse tomographic reconstruction. In accordance with this method, projections are measured from transmission, emission and/or reflection data collected from a subject under study. A system matrix is then constructed to relate the measured projections to a reconstructed tomographic image, this system matrix comprising vector-columns respectively representative of measured projections. The system matrix is factored by singular value decomposition, wherein this factored system matrix includes a diagonal matrix of singular values. The factored system matrix is pseudo-inverted to produce a pseudo-inverted system matrix in which the singular values are presented as a non-increasing sequence called singular value spectrum.
The pseudo-inverted system matrix is regularized by truncating the singular value spectrum at an index T and removing small singular values /si, i = T +1, ..., M, M being the total number of image voxels, or modifying the singular value spectrum to diminish the effect of very small singular values. Finally, a new reconstructed image is evaluated by summing a previous image estimate and one vector-column of the pseudo-inverted regularized system matrix corresponding to a last event recorded. In operation, reconstruction of the tomographic image involves summation of the previous image estimate and a current update from a new measured event or additional single projection data.
The method of the invention presents, amongst others, the following advantages:
- it allows real-time reconstruction of projection data;
- it allows event-by-event updating of a tomographically reconstructed image;
- it allows individual projection data to be reconstructed independently to update an existing image;
- it allows real-time updating and display of an image;
- it allows to monitor tomographic data acquisition in real-time;
- it avoids storage of measured projection data (either in list mode, histogram or rebinned sinogram); and - it avoids data rebinning into projection vectors, since the exact geometry of the measuring instrument can be utilized.
The above as well as other advantages and features of the present invention will become more apparent upon reading of the following non restrictive description of a preferred embodiment thereof, given as illustrative example only with reference to the accompanying drawings.
BRIEF DESCRIPTION OF THE DRAWINGS
In the appended drawings:
Figure 1 is an illustration of the various steps involved in the tomographic image reconstruction method by filtered backprojection (FBP);
Figure 2 is a comparison of the same PET tomographic data reconstructed by FBP (left) and a statistical iterative method maximizing likelihood with an expectation maximization algorithm (right), in which the image represents a radiotracer uptake in two 25 g mice, 3 minutes after injection of NCi 18F-FDG, and in which a 20 seconds scan 0105 counts) was obtained with the Sherbrooke animal PET scanner;
Figure 3 is a schematic diagram of dedicated hardware for performing real-time tomographic image reconstruction based on the matrix pseudo-inverse method;
Figure 4 is a singular value spectrum of the system matrix for an animal PET scanner, 64x64 pixel image;
Figure 5 is a sequence of phantom images demonstrating the use of real time TSVD reconstruction, wherein estimated radioactivity distribution progress is shown with one plane (64x64 pixel images), approximately 5400 additional counts are accumulated with each next image estimate shown, and less than 1/3 of the singular value spectrum is utilized; and Figure 6 is a s equence of phantom images demonstrating the use of real time TSVD reconstruction, wherein estimated radioactivity distribution progress is shown with one plane (64x64 pixel images), projections at two additional angles are included into each next image estimate shown, and less than 1/3 of the singular value spectrum is utilized.

DETAILED DESCRIPTION OF THE PREFERRED EMBODIMENT
A different approach that was discussed in the literature [21 ] some time ago, but has not yet gained popularity for tomographic reconstruction is the matrix pseudo-inverse method, whereby the system matrix relating the measured projections to the reconstructed tomographic image must be inverted.
Its use has been hindered by the inverse problem ill-conditioning and the computational burden resulting from the size of the matrix. Matrix inversion can be performed based on singular value decomposition (SVD) of the system matrix with some regularization to diminish the effect of small singular values. The application to the tomographic reconstruction problem is described below, followed by examples of practical hardware implementation of the technique for real-time image reconstruction.
Matrix Pseudo-Inverse Tomographic Reconstruction The image reconstruction problem in tomography can be written in the matrix form Px = b (4) where N
P = p~ : i =1,..., N; ~ ply =1, j =1,..., M (5) t=I
is the system matrix, b = {b1: i =1,..., N~ is the vector-column of measured projections (sinogram), x = {xk : k =1,...,M} is the vector-column (image) that has to be found, N and M are the total number of tubes-of response and image pixels (or voxels), respectively. For example, the matrix P may be obtained in PET for a given scanner geometry with the approach described in [18].
5 Any matrix P can be factored [22] into P = UWV T (g) where U = ~u,~ : i =1,...,N; j =1,...,M~ and V = ~v~~ : i, j =1,...,M~ are orthogonal 10 matrices, W = ~,u~~ : i, j =1,...,M; elf = 0 if i ~ j~ is a diagonal matrix containing singular values ,u; ---,u;;, i =1,...,M . Factored matrix representation (6) is referred to as singular value decomposition (SVD).
One may find a least-squares solution of (4) using x = P+b (7) where P+ =VW+UT (8) is the pseudo-inverse of P [22]. Singular values are usually presented as a non-increasing sequence called singular value spectrum. Some singular values can be very small (or even zeros if the matrix is singular).
The condition number max fc;
c - ' (9) min,u;
t of the system matrix can be very high (even infinity for a singular matrix).
Thus, the inverse problem of tomographic image reconstruction is ill-conditioned and the solution exploiting pseudo-inverse (8) will be very sensitive to noise.
One way of solution regularization is to truncate the singular value spectrum at some index T and remove small singular values ,u~, i = T +1,...,M
from the solution expansion, leaving T N
x = xk : xk = ~vki,ua 1 ~u jtb j; k =1,...,M (10) i=1 j=1 which is referred to as the truncated SVD (TSVD) solution.
Another regularization approach would be to modify the singular spectrum in order to diminish the effect of the inversion of very small singular values (and zeros if any):
x = zk : .zk = ~vka f ~,u~ ~ 1 ~ a jab j; k =1,...,M (11 ) i=1 j=1 where f (,u~ ~ is the modified singular value. Such solution is referred to as the modified SVD (MSVD) solution. An example of MSVD can be found in [23].

Real-Time Reconstruction In this section let us concentrate on the TSVD solution. Rearranging the summation order in (10) we get T N N T _ xk = ~ ~ vki~i lujibj = ~ ~vki~i luji bj, k=I,...,M (12) i=lj=I j=1 i=I
Let P+ be the "truncated pseudo-inverse":
_ _ T
p+ _ ~Pkj v Pkj = ~vkil~i lujie k =1,...,M; j =1,...,N~ (13) i=I
Then the TSDV solution in the matrix form is z = P+b (14) Let b(t~ be a sinogram containing a total of t counts:
N
b(t~= bi(t~---bi : t=1,...,N; ~bi =t (15) i=1 and x(t~ _ ~xk (t~ : k =1,..., M~ is the respective solution given for b(t~ by (14). Let Ob(S) _ Obi : i =1,...,N~ is such that O,i=1,...,s-1 Ob;--_- l,i=s (16) 0, i = s + 1,..., N
and b(S)(t) ---- b(t -1~+ 0 b(S) (17) In other words, b(S)(t)= ~b~s)(t~: i =1,...,N~ is a sinogram, which differs from b(t-1) by only one count in a bin with index s:
b; (t-1~, i =1,...,s-1 b~s)(t)= bi(t-1~+1, i =s (18) b; (t-1~, i =s+1,...,N
The solution given by x~s)(t~= p+b~s)(t~ (1 may be found as follows. We have xk(t~= ~Pkjbjs)(t~= ~Pkjbi(t-l~+Pks~bs(t-1~+l~+
j=t j=t (20) + ~ Pkjbt(t-1~= ~Pkjb~(t-1~+Pks~k=1~...,M
j=s+1 j=I

Thus, the new image estimate is simply a sum of the previous image estimate and one column (vector-column) of matrix P+:
x(S)(t~= ~ks)~t~: xks)(t)= xk(t-1~+ per, k =1,...,M} (21) It is natural to start from the blank image x = ~xk : k =1,...,M; xk = 0} and use (21 ) afterwards to obtain the current image estimate in real-time based on the previous image and a column of P+ defined by the sinogram bin index where the last event was registered.
MSVD solution can be obtained in real time as well using (21 ) and taking into account that "modified pseudo-inverse" P+ is given in this case by:
M _ P+ _ ~pk~ ~ pk~ --- ~ vkaf (l~~ ) 1 u~j ; k =1,..., M; j =1,..., N (22) a=~
The proposed approach is general and can be applied with any number of projection bins and image pixels (or voxels), for 2D as well as for 3D-reconstruction geometry.
Corrections of projection data Physical and instrumental determinants that affect the measurement of projection data must be incorporated into the image reconstruction. Such factors as the spatially variant system response and model of the emission and detection processes can be included into the computation of the system matrix.

Other factors related to the specific instrument being used and the subject under study, such as the normalization for detector efficiency or the correction for signal attenuation in tissue (e.g., photon attenuation in the case of emission tomography), can be included by multiplying the pseudo-inverse matrix elements 5 with the proper factors.
The updated equation (21 ) will change slightly in the case of correction on the fly:
10 x(S)(t~=~cks)~t~: xks)~t~=xk~t-1~+FS x per, k=1,...,M} (23) where FS is the correction factor assigned to the tube-of-response s.
Similarly, random coincidence events in PET could be corrected for by 15 subtracting one column of matrix P+ from the previous image estimate:
x(S)~t~=~cks)~t~: xks)~t~=xk(t-1~-P~,k=1,...,M} (24) or, instead, by skipping the next column addition given by (21 ) when the next coincident event is registered. In the case of correction on the fly, (24) would be replaced by:
x(S)(t)=~ks)(t): xks)(t)=xk~t-1~-FS x per, k=1,...,M} (25) Hardware implementation The reconstruction of tomographic images using the proposed method has been shown to be amenable to a simple sum of two vectors representing the previous image estimate and the current update from a new measured event or additional single projection data. The latter vector may be multiplied by a normalization factor accounting for physical or instrumental corrections.
Computation must be performed in floating point arithmetic to satisfy the required precision.
For real-time reconstruction, the computation is subject to the additional constraint that the duration of the entire calculation must be less than the average time between successive events. This includes fetching the column of P+defined by the bin index of the current event, multiplying by a correction factor if necessary, summing the M elements of the resulting vector to the previous image x = ~xk : k =1,...,M~, and storing the updated image.
Two different acquisition modes must be distinguished here, depending on the way projection data are being generated by the tomographic scanner. In the first mode, the measuring instrument sequentially supplies projection data, which are measured only once in each projection bin. X-ray CT scanners and ultrasound probes fall into this category. In such cases, the average time between events is determined by the scanning speed of the instrument, and is typically in the millisecond range or more. Computation in such case can be performed using general-purpose computers.
In the second mode of acquisition, events are registered randomly in all available projection bins. Instruments operating in counting mode, such as SPECT and PET scanners, fall into this category. Assuming a peak event rate of 106 counts/sec, the entire computation would have to be accomplished in less than one microsecond in such case. Even for a small size 2D image (e.g., 64x64 pixels, M=4096), such computational speed can hardly be reached with the current general-purpose computer technology or even with specialized digital signal processors (DSPs), unless massively parallel processing is employed.
The system matrix P+, having dimensions N x M, must also be loaded into memory for fast access. Therefore, extending the reconstruction to larger image sizes or 3D imaging geometry makes such implementation impractical with the current computer technology.
Dedicated hardware with highly parallel processors and high storage capacity must be designed to efficiently perform the required computation. The principle of a real-time processing unit to perform tomographic image reconstruction based on the matrix pseudo-inverse method is displayed schematically in Figure 3. Such hardware can easily be implemented using programmable logic (e.g., Field Programmable Gate Arrays or FPGA) or application specific integrated circuits (ASIC) by persons knowledgeable in the field.
Results As an example of the proposed reconstruction method, the singular value spectrum of the system matrix P for the above mentioned animal PET scanner [24] and image of 64x64 pixels is shown in Figure 4. The matrix condition number is c = 4441,5 in this case.
Sequence of phantom images demonstrating the use of real time TSVD
reconstruction based on (18) is presented in Figures 5 and 6. Data acquired with a phantom of 110 mm diameter made of Lucite and having holes of diameter 2, 3.4, 6.7, 9.7, 13, 15.8, 20.3, 22.7 mm located on a circumference at a distance of 28 mm from the center and filled with 18F-fluorodeoxyglucose were reconstructed. In Figure 5, every image has approximately 5400 additional counts as compared to the previous image estimate. Figure 6 shows image estimates obtained as projections at consecutive angles. In both cases, less than 1/3 of the singular value spectrum was utilized.
Discussion of application in medical ima in SVD of the system matrix, apart from precise numerical diagnostics of the tomographic reconstruction ill-conditioning with a given detection system geometry, provides a linear and very fast reconstruction means.
Singular value spectrum truncation allows separation of the signal and noise at the reconstruction step. Index T sets the trade-off between noise and resolution. Truncation of the singular spectrum is not the only way of solution regularization. Spectrum modification without truncation may be an appropriate regularizing approach in some situations. It is possible as well that in some special situations when the system matrix corresponding to a particular tomographic system is sufficiently well conditioned no SVD pseudo-inversion is necessary and direct matrix inversion is feasible.
TSVD reconstruction has some drawbacks: negative values in the image estimate, streak artifact with low-count images if index T is lower than a certain threshold value or noise artifact if T is higher than the threshold value. But these features (except for the noise artifact with high T, which is analogous to the artifact developing with high iteration numbers when unconstrained iterative image estimation is performed) are also common to the filtered backprojection image reconstruction, the most popular method used in medical practice today.
TSVD (as well as MSVD) also has very attractive benefits:
1 ) spatially variant system response and model of the signal emission and detection processes can be easily included into image reconstruction through the system matrix;
2) data rebinning is not necessary since the geometry of a given system is utilized;
3) the resolution in reconstructed images may be set based on the spatial resolution analysis in reconstructed images by varying truncation index T (or modifying singular value spectrum in MSVD case);
4) noise amplification may be controlled by varying truncation index T (or modifying singular value spectrum in MSVD case);
5) image reconstruction can be performed in real time, on an event-by-event or projection-by-projection basis, allowing for instant visualization of the radioactivity distribution while the subject is being scanned; and 6) since measured projection data can be reconstructed "on the fly" as soon as they are acquired by the scanner, storage of the measured data or intermediary calculation results is unnecessary, which allows for optimum compression of the tomographic data.
Although the present invention has been described hereinabove by way of a preferred embodiment thereof, this embodiment can be modified at will, within the scope of the present invention, without departing from the spirit and nature of the subject invention.
RAfArr~nrc~ c-[1] J. Radon, "On the determination of functions from their integrals values along certain manifolds", IEEE Trans. Med. Imag., vol. MI-5, no. 4, pp. 170-176, 1986.
(Translated from the original text that appeared in Berichte Saechsische Akademie des Wissenschaften, vol. 29, pp. 262-277, 1917) [2] G.N. Hounsfield, "Method and apparatus for measuring X- or o-radiation absorption or transmission at plural angles and analyzing the data", U.S.
Patent 3,778,614, 1973.
5 [3] G.N. Hounsfield, "Computerized transverse axial scanning (tomography).
1.
Description of system", Brit. J. Radiol., vol. 46(552), pp. 1016-22, 1973.
[4] A.M. Cormack, "Representation of a function by its line integrals, with some radiological applications", J. Appl. Phys., vol. 34, pp. 2722-2727, 1963.
[5] A.M. Cormack, "Representation of a function by its line integrals, with some 10 radiological applications. II", J. Appl. Phys., vol. 35, pp. 2908-2913, 1964.
[6] A.M. Cormack, "Reconstruction of densities from their projections, with applications in radiological physics", Phys. Med. Biol., vol. 18(2), pp. 195-207, 1973.
[7] L.A. Shepp, B.F. Logan, 'The Fourier reconstruction of a head section", IEEE
15 Trans. Nucl. Sci., vol. NS-21, pp.21-43, 1974.

[8] G.T. Herman, Image Reconstruction from Projections, Academic Press, New York 1980.

[9] S.R. Dean, The Radon Transform and Some of Its Applications, John Wiley, New York, 1983.
20 [10] F. Natterer, The Mathematics of Computerized Tomography, Wiley, New York, 1986.
[11] A.C. Kak, M.Slaney, Principles of Computerized Tomographic Imaging, IEEE Press, New York, 1988.
[12] J.A. Fessler, J.M. Ollinger, Signal Processing Pitfalls in Positron Emission Tomography, Technical Report 302, Comm. and Sign. Proc. Lab., Dept. of EECS, Univ. of Michigan, 48109-2122, Ann Arbor, MI, 1996.
(http://www.eecs.umich.edu/~fessler).
[13] Y. Vardi, L.A. Shepp, L. Kaufman, "A statistical model for positron emission tomography", J. Am. Stat. Assoc., vol. 80, pp. 8-37, 1985.
[14] J.M. Ollinger and J.A. Fessler, "Positron-Emission Tomography", IEEE
Signal Processing Mag., vol. 14(1 ), pp. 43-55, 1997.

[15] G.T. Herman, L.B. Meyer, "Algebraic reconstruction techniques can be made computationally efficient", IEEE. Trans. Med. Imag., vol. 12, no. 3, pp.
600-609, 1993.
[16] X.-L. Xu, J.-S. Liow, S.C. Strother, "Iterative algebraic reconstruction algorithms for emission computed tomography: A unified framework and its application to positron emission tomography", Med. Phys., vol. 20(6), pp. 1675-1684, 1993.
[17] A.J. Rockmore, A. Macovski, "A maximum likelihood approach to emission image reconstruction from projections", IEEE Trans. Nucl. Sci., vol. 23, pp.

1432, 1976.
[18] V.V. Selivanov, Y. Picard, J. Cadorette, S. Rodrigue, R. Lecomte, "Detector response models for statistical iterative image reconstruction in high resolution PET", IEEE Trans. Nucl. Sci., vol. 47, 2000 (to appear) [19] V. Selivanov, D. Lapointe, M. Bentourkia, R. Lecomte, "Cross-validation stopping rule for ML-EM reconstruction of dynamic PET series: Effect on image quality and quantitative accuracy", 1999 IEEE Nuclear Science Symposium &
Medical Imaging Conference Record, in press. (CD, MIC Sessions, M10-18) [20] E. Levitan, G.T. Herman, "A maximum a posteriori probability expectation maximization algorithm for image reconstruction in emission tomography", IEEE
Trans. Med. Imag., vol. 6, no. 3, pp. 185-192, 1987.
[21] Y.S. Shim, Z.H. Cho, "SVD pseudoinversion image reconstruction", IEEE
Trans. ASSP, vo1.29, no.4, pp.904-909, 1981.
[22] G. Strang, Linear algebra and its applications, Harcourt Brace Jovanovich, 1980, pp. 137-145 [23] M.F. Smith, "Generalized matrix inverse reconstruction for SPECT using a weighted singular value spectrum", IEEE Trans. Nucl. Sci., vol. 43, no. 3, pp.
2008-2017, 1996.
[24] R. Lecomte, J. Cadorette, P. Richard, S. Rodrigue, D. Rouleau, "Design and engineering aspects of a high resolution positron tomograph for small animal imaging", IEEE Trans. Nucl. Sci., vol. 41, pp. 1446-1452, 1994.

Claims

The embodiments of the invention in which an exclusive property or privilege is claimed are defined as follows:

1. A method for conducting real-time matrix pseudo-inverse tomographic reconstruction, comprising:
measuring projections from transmission, emission and/or reflection data collected from a subject under study;
constructing a system matrix relating the measured projections to a reconstructed tomographic image, said system matrix comprising vecto-columns respectively representative of measured projections;
factoring by singular value decomposition the system matrix, wherein the factored system matrix includes a diagonal matrix of singular values;
pseudo-inverting the factored system matrix to produce a pseudo-inverted system matrix in which the singular values are presented as a non-increasing sequence called singular value spectrum;
regularizing the pseudo-inverted system matrix by truncating the singular value spectrum at an index T and removing small singular values µi, i = T +1, ..., M, M being the total number of image voxels, or modifying the singular value spectrum to diminish the effect of very small singular values;
and estimating a new reconstructed image by summing a previous image estimate and one vector-column of the pseudo-inverted regularized system matrix corresponding to a last event recorded;
whereby reconstruction of the tomographic image involves summation of the previous image estimate and a current update from a new measured event or additional single projection data.