CA2360656A1 - Real-time image reconstruction for computed tomography systems - Google Patents

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 system 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, digital signal processors, dedicated programmable logic devices or application specific circuits. 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 pseudo-inverted and regularized 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 optimal 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 and apparatus for conducting matrix pseudo-inverse tomographic reconstruction, in particular but not exclusively in real-time.

2. Brief description of the current technology:
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 Radon in 1917, hence the denomination inverse-Radon transform. However, it was only in the early 1970's, with the invention of the X-ray computed tomographic (CT) scanner and the development of reconstruction algorithms, that tomographic imaging came into widespread use. 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
(three-dimensional) ultrasound (US), were the result of the application of tomographic principles.
A. Review of current tomo4raphic reconstruction algorithms Many tomographic reconstruction algorithms are disclosed in the art and a comprehensive description of these can be found in the literature. 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 improvement of image estimates at every iteration.
A.1 Filtered Back Pro~iection The most common tomographic reconstruction method nowadays is the filtered back projection (FBP) algorithm, which makes use of the projection-slice theorem, see generally A.C. Kak & M. Slaney, Principles of Computerized Tomographic Imaging, IEEE Press, New York, 1988. The algorithm generally 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 back projecting the result onto the image grid:
z =Back-Project{F~'~c(f)xF,~b'~~ (1) where b' is a j-th projection vector (j=1,...,A), A is a total number of projection angles, F, 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 in the simplest noise-free case, 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 filtering step, which is commonly pertormed as a multiplication with the ramp filter in Fourrier space, can also be carried out as a convolution of the projection vector in the spatial domain with a filter kernel that is the inverse Fourier transform of the ramp filter.
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 scanner geometry), or must be rebinned with equal spacing to accommodate the convolution operation (commonly 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 artifacts are present; and - event-by-event reconstruction involves some elaborate computation to align, rotate and scale the single event filters before summation to the image matrix.
A.2 Algebraic Reconstruction In the art there are known other tomographic image reconstruction methods called algebraic reconstruction techniques. They aim at solving the system of linear equations:
Pl lxl + P12x2 '~' ... + plMxM = b1 P21x1 +P22x2 +...+ p2MxM = b2 (2) PNlx1 + PN2x2 +... + pNMxM = bN
where bt : i =1,..., N are the measured projections (sinogram), x = {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 preferably (though not always) converge to the minimum norm least-squares solution as the iteration number goes to infinity, see in particular 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. Some of the drawbacks of the iterative algebraic methods are:

5 - they are very computer intensive and have large memory requirements;
- given the inherently iterative nature of the algorithms, real-time reconstruction is virtually not implementable;
- the need to optimize the iteration number and the order in which projections are being utilized in the reconstruction;
- the non-linearity of the solution due to the iterative nature of reconstruction algorithm, which results in complex quantitative dependence of the image estimate on input data; 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 account of the stochastic nature of the measured data in emission tomography and the physical modeling of the detection process. In general, statistical iterative reconstruction methods yield images with a greatly superior visual quality when compared to FBP reconstructed images, especially in the case of low projection statistics. The absence of negative values and proper modeling of the detector response helps to avoid several of the artifacts present in FBP
images. Unfortunately, the superior image quality does not necessarily translate into better quantitation accuracy. The use of statistical reconstruction methods also raises other practical problems that still hamper their utilization outside the research environment. In addition to most of the drawbacks that are common with iterative algebraic methods, statistical methods also give rise to the following problems:
- in the case of unconstrained iterative estimation, the need to optimize the number of iterations: If the number of iterations is too low, the spatial resolution is sub-optimal; if the number of iterations is too high, the spatial resolution is artificially enhanced at the expense of unacceptable noise amplification;
- in the case of so called Bayesian approach the necessity to optimize an extra parameter, called prior, that controls noise or image smoothness (see E.
Levitan & G.T. Herman, "A maximum a posteriors probability expectation maximization algorithm for image reconstruction in emission tomography", IEEE Trans. Med. Imag., vol. 6, no. 3, pp. 185-192, 1987); and - the non-linearity of the solution due to the iterative nature of the reconstruction algorithm, which results in a complex quantitative dependence of the image estimate on input data.
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 reconstructing an image, formed of an array of pixels (or voxels), from a series of measured individual tomographic projection data, comprising creating a system matrix defining vectors relating the measured individual projection data to the pixels (or voxels) of the image, associating one of the vectors of the system matrix to each measured individual projection data and, for each measured individual projection data of the series, updating the image by adding to a previous image estimate the vector of the system matrix associated to the measured individual projection data.
According to preferred embodiments of the image recontructing method:
- the image reconstructing method comprises conducting the vector associating and the image updating in real-time while the individual projection data are measured;
- creating a system matrix comprises factoring a system matrix defining vectors relating the measured individual projection data to the pixels of the image, and inverting the factored system matrix;
- the factored system matrix presents the form P=UWVT , and inverting the system matrix comprises pseudo-inverting the factored system matrix of the form P=UWVT to obtain a pseudo-inverse system matrix of the form P+=VIiV+UT where U and V are orthogonal matrices, W is a diagonal matrix containing singular values, W'' is the inverted diagonal matrix hV, and T
denotes transpose of matrices V and U;
- the image reconstructing method comprises truncating the pseudo-inverse system matrix to remove small singular values;
- the image reconstructing method comprises modifying a singular value spectrum to reduce the effect of inversion of very small singular values, wherein modifying the singular value spectrum comprises replacing the singular values in the diagonal matrix W''' by a function of these singular values;
- the measured individual projection data comprise events detected through tubes-of-response of a tomographic scanner, the system matrix comprises a number of columns equal to the number of tubes-of-response and a number of rows equal to the number of pixels of the image, the columns of the system matrix define respective vector-columns relating the events detected through the tubes-of response to the pixels of the image, each tube-s of response is associated to a respective vector-column of the system matrix, associating one of the vectors of the system matrix to each measured individual projection data comprises associating each detected event to the vector-column of the system matrix associated to the tube-of response through which the event has been detected, and updating the image comprises, for each detected event, adding to the previous image estimate the vector-column of the system matrix associated to the tube-of response through which the event has been detected;
- the image reconstructing method comprises, prior to adding to the previous image estimate the vector-column, scaling that vector-column with a correction factor; and - the image reconstructing method comprises obtaining the previous image estimate by adding together all vectors of the system matrix associated to all measured individual projection data prior to measurement of the individual projection data under consideration.
The present invention also relates to an apparatus for reconstructing an image, formed of an array of pixels, from a series of measured individual tomographic projection data and a system matrix defining vectors relating the measured individual projection data to the pixels of the image. This image reconstructing apparatus comprises means for associating one of the vectors of the system matrix to each measured individual projection data, and means for updating the image for each measured individual projection data of the series.
The latter updating means comprises an adder for adding to a previous image estimate the vector of the system matrix associated to the measured individual projection data.

Advantageously, this image reconstructing apparatus includes means for operating the associating means and the updating means in real time while the individual projection data of the series are measured.
Still further in accordance with the present invention, there is provided an apparatus for reconstructing an image, formed of an array of pixels, from a series of measured individual tomographic projection data and a system matrix defining vectors relating the measured individual projection data to the pixels of the image. This image reconstructing apparatus comprises a generator of a projection data index in response to each measured individual projection data, a memory unit storing the system matrix and comprising an address input supplied with the projection data index from the generator and a data output delivering a vector of the system matrix corresponding to the projection data index and, for updating the image for each individual projection data, an adder for adding to a previous image estimate the vector of the system matrix delivered on the data output of the memory unit.
According to preferred embodiments of the image reconstructing apparatus:
- the image recontructing apparatus comprises means for operating the generator, memory unit and adder in real time while the individual projection data are measured;
- the system matrix is a pseudo-inverse system matrix of the form P+=VW+UT, which results from a peuso inversion of a system matrix having the form P=UWVT, where U and V are orthogonal matrices and W is a diagonal matrix containing singular values, W+ is the inverted diagonal matrix W, and T denotes transpose of matrices V and U;

- the image reconstructing apparatus comprises means for truncating the pseudo-inverse system matrix to remove small singular values;
- the image reconstructing apparatus comprises means for modifying a singular value spectrum to reduce the effect of inversion of very small 5 singular values, wherein the means for modifying the singular value spectrum comprises means for replacing the singular values of the diagonal matrix INS
by a function of these singular values;
- the measured individual projection data comprise events detected through tubes-of-response of a tomographic scanner, the memory unit comprises a 10 number of columns of memory locations equal to the number of tubes-of response and a number of rows of memory locations equal to the number of pixels of the image, the columns of the memory unit store respective vector-columns relating the events detected through the tubes-of response to the pixels of the image, each tube-of-response is associated to a respective column of the memory unit, and the adder has a first input receiving the previous image estimate and a second input receiving the vector-column stored in the column of the memory unit associated to the tube-of response through which the event has been detected;
- the image reconstructing apparatus further comprises a multiplier which, prior to adding to the previous image estimate the vector-column, multiplies this vector-column with a correction factor; and - the image reconstructing apparatus comprises means for obtaining the previous image estimate by adding together all vectors of the system matrix associated to all measured individual projection data prior to measurement of the individual projection data under consideration.
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 data rebinning into projection vectors, since the exact geometry of the measuring instrument can be utilized; and - it may avoid storage of measured projection data (either in list mode, histogram or rebinned sinogram) if further processing of the measured projection data is unnecessary.
The foregoing and other objects, 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 a is a schematic diagram of the detection unit of a PET scanner;
Figure 1 b is a schematic diagram of the detection unit of a PET scanner showing, for a fan-beam response, tubes-of response between one detector paired with a plurality of other detectors;
Figure 1c is a schematic diagram of the detection unit of a PET scanner showing, for linear sampling, parallel tubes-of-response between pairs of detectors in a detector array;
Figure 2a is a schematic diagram of a phantom used for generating test images in a PET sanner;
Figure 2b is an illustration of the various steps involved in a real-time tomographic image reconstruction using the matrix pseudo-inverse method;
Figure 3 is a schematic diagram of a 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 and a 64x64 pixel image;
Figure 5 is a first sequence of phantom images demonstrating the use of real-time TSVD reconstruction; and Figure 6 is a second sequence of phantom images demonstrating the use of real time TSVD reconstruction with projections acquired at eight (8) consecutive angles.

DETAILED DESCRIPTION OF THE PREFERRED EMBODIMENT
An approach that was disclosed some time ago (see Y.S. Shim & Z.H.
Cho, "SVD pseudoinversion image reconstruction", IEEE Trans. ASSP, vo1.29, no.4, pp.904-909, 1981), but has not yet gained popularity for tomographic reconstruction is the matrix pseudo-inverse method, whereby the system matrix relating the tomographic image grid (pixels) to the measured projections must be pseudo-inverted. Its use has been suggested in situations where only a limited number of projections are required. To provide complete number of projections, which is necessary in any practical application of tomographic medical imaging, its use has been hindered by the computational burden resulting from the size of the matrix and ill-conditioning of the system matrix which has to be accounted for in a manner preferably independent of the image being reconstructed. Matrix pseudo-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. Its application to the tomographic reconstruction problem is described below.
Matrix Pseudo-Inverse Tomograahic Reconstruction The image reconstruction problem in tomography can be written in the matrix form Px = b (3) where N
P = p J : i =1,..., N; ~ p~~ =1, j =1,..., M (4) t=i is the system matrix, b = {b~ : i =1,..., N~ is a vector-column of measured projections (sinogram), x = {xk : k =1,...,M} is a vector-column (image) that has to be found, and N and M are the total number of tubes-of response and image pixels (or voxels), respectively. As known to those of ordinary skill in the art, the matrix P may be obtained in PET for a given scanner geometry with the approach described in 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, No. 3, pp. 1168-1175, 2000.
Any matrix P can be factored into P = UWV T (5) where U = {ul~ : i =1,..., N; j =1,..., M and V = {v;~ : i, j =1,..., M are orthogonal matrices, W = f ,ulf : i, j =1,..., M; ~Z~ = 0 if i ~ j is a diagonal matrix containing singular values ,u1 ---- ,u11, i =1,..., M . T denotes matrix transpose.
Factored matrix representation (5) is referred to as Singular Value Decomposition (SVD). An algorithm to compute a matrix SVD is disclosed in W. H. Press & B. P.
Flannery, S. A. Teukolsky, W. T. Vetterling, "Numerical recipes in C: The art of scientific computing", Cambridge University Press, 1988, pp. 60-72.
In the case where measured data are noisy and/or Equation (3) is overdetermined/underdetermined, and the exact solution of Equation (3) does not exist, one may find a minimum norm least-squares solution of Equation (3) using z = P+b (6) where P+ = VW +UT (7) 10 is the pseudo-inverse of P. tN+ is the pseudo-inverse diagonal matrix hV.
Singular values are usually presented as a non-increasing sequence called singular value spectrum. Some singular values can be very small (or in some cases equal to zero if the matrix is singular).
15 The condition number max,u;
c = ' (8) min ~1 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 any solution attempting to exploit the pseudo-inverse of matrix P (7) will be very sensitive to noise.
One way of regularizing the solution 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 _1 N
x = zk : zk = ~vktf~~ ~u~tbJ ~ k =1,...,M (9) t=t >=1 which is referred to as the truncated SVD (TSVD) solution.
Another regularization approach would be to modify the singular value spectrum in order to diminish the effect of the inversion of very small singular values (and zeros if any):
M N
z = zk : xk = ~ vk; f (,tt; y u;,b~; k =1,..., M (10) i_t where f (,u~ ~ is a function of the singular value. Such a solution is referred to as the modified SVD (MSVD) solution. An example of the MSVD can be found in R. Lecomte, J. M.F. Smith, "Generalized matrix inverse reconstruction for SPECT using a uveighted singular value spectrum", IEEE Trans. Nucl. Sci., vol.
43, no. 3, pp. 2008-2017, 1996.
Real-Time Reconstruction Rearranging the summation order in Equation (9) we get:

xk = ~ E vk~W lu~tb~ _ ~ ~vkr~t lu~~ b~; k=1,...,M (11) f=Ij=1 ~=1 i=~
Let P+ be the "truncated pseudo-inverse":
_ T
p+ - Pk~ ~ Pk~ _ ~vkil~i ~u~t~ k =1,...,M; j =1,...,N (12) t=t Then the TSDV solution in the matrix form is:
z = P+b (13) Let b(t) be a sinogram containing a total of t counts:
N
b(t) = b~ (t) --__ b; : i =1,..., N; ~ b1 = t ( 14) c=i and x(t) _ ~xk (t~ : k =1,..., M is the respective solution given for b(t) by Equation (13). Let Ob(S) _ {O b; : i =1,...,N such that 0, i =1,..., s -1 Obl-~-= l,i=s (15) O,i=s+1,...,N
and b(s) (t~ --- b(t -1~+ 0 b(S) (16) 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, the count being located in a bin with index s:
bZ(t-1~, i =1,...,s-1 b~s)(t)= b1 (t-1~+l,i=s (17) b; (t-l~,i=s+1,...,N
Assuming that x~s~(r) stands for the image estimate obtained using input data b~'~(t~, then the solution given by z(S)(t~= P+b(S)(t) (18) may be found as follows. Given that the sinogram b~s~(t~ differs from b(t-1~
by only one (1 ) count, we have:
N s-1 xx(t~- ~pkjbjs~~(t~- ~pkjbj(t-1)+ pksfbs(t-1)+lJ
j=1 l=1 N _ N (19) + ~Pkjbj(t-1~=~Pkb;(t-l~+P~~ k=1~...,M
j=s+1 j=1 Thus, the new image estimate is the sum of the previous image estimate and one column (vector-column) of matrix P+
x(S)(t)= ~xks)(t): xks)(t)= xk(t-1~+ per, k =1,...,M (20) The natural starting point for the described image reconstruction is the blank image x(0~= ~xk (0~: k =1,..., M; xk (0~= 0 , using Equation (20) subsequently to obtain a current image estimate in real-time based on the previous image and the column of P+ defined by the sinogram bin index where the last event was registered.
The MSVD solution can be obtained in real time as well using Equation (20) and the modified pseudo-inverse matrix P+, given in this case by:
_ M
- hkj ~ Pkj = ~ vki.~(f~i ~ ji ~ k =1,..., M; j =1,..., N (21 ) i=1 where ~ is a function of other than ~ _ ~' ~' t ~ T
f (,u; ,u; ( f (~; , which was used 0, f > T
for deriving the TSVD).
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, with scanners collecting complete or incomplete sets of projections, for non-singular or singular system matrices. The proposed approach is not limited to PET, but can also be applied with any tomographic imaging system, such as SPECT, CT, etc., as long as a system matrix P relating the image pixels to the measured projection data can be obtained and inverted.
In the following preferred embodiment a PET scanner having 2D reconstruction geometry is disclosed.

Referring now to Figure 1a, a radioisotope which decays by positron (~i+) emission is injected into (or inhaled by) a patient (not shown). The spatial and temporal distribution of the radioisotope within the patient will then depend on the manner in which the organ or tissue being scanned manages the 5 radioisotope both biochemically and physiologically. Once within the patient, the positrons escape from the radioisotope and travel a short distance, colliding with electrons of nearby atoms. The collision of a positron and electron results in an annihilation reaction 2 whereby the mass of the positron and electron are converted into a pair of coincident 511 keV photons 4 which are emitted at 180°
10 to one another (note that a very small portion of the annihilation reactions result in the emission of three (3) annihilation photons but these are statistically unimportant). The annihilation photons easily escape the patient's body and can be recorded by detectors 6 arranged in a circular array, or any other convenient geometry, around the patient.
The detectors 6 which are used in a PET scanner such as the one generally designated by the reference 1 in Figure 1 a, typically consist of scintillation crystals 8 of bismuth germanate (BGO) which convert the gamma rays into light photons. The light photons are converted in the detectors 6 (for example, photomultiplier tubes or avalanche photodiodes) into electrical signals 10 which can then be transferred to the PET scanner's processing subsystem 12. The detectors 6 are typically arranged in rings such as14, generally with a number of rings such as 14 distributed along an axis 16 making up the PET
scanner.
As stated above, the coincident rays 4 are emitted simultaneously and at 180° to one another and, therefore, by detecting the occurrence of simultaneous (i.e. within 5-15 nanoseconds typically) reception of coincident rays 4 at two different detectors 6 (known as a coincident event), the annihilation reaction, and therefore the radioisotope, can be localised within the patient along a line drawn between the two detectors 6 having sensed the two coincident rays 4.
As higher concentrations of a positron emitting radioisotope give rise to a larger number of annihilation reactions 2 and therefore coincident events, the processing subsystem 12 can process the coincident events and reconstruct an image of the concentration and location of the radioisotope within the patient.
Referring now to Figure 2a, the process of real-time tomographic image reconstruction in its preferred PET embodiment will be further described. To model the effects of a positron emitting substance injected into a patient while at the same time providing an image which is easily recognizable, a cylindrical phantom 18 is used. The phantom 18 is preferably constructed of Lucite~
(plexiglass material) or an equivalent material with a series of holes 20 of different diameters for accepting a positron emitting substance 22, for example 18F-fluorodeoxyglucose (FDG). The phantom 18 is placed within the PET
scanner (not shown) for recording coincident annihilation photon emissions. At time TO the positron emitting substance 22 emits positrons which collide with electrons of nearby atoms (neither shown) resulting in an annihilation reaction, the generation of two (2) annihilation photons 24180° apart from each other and the detection of a new coincident event 26 by a pair of detectors 6.
Referring back to Figure 1a, as stated above, the BGO detectors 6 of the PET scanner 1 are typically arranged in rings 14. The region within the center of the PET scanner 1 between a given pair of detectors 6 is defined as a tube-of response 28. Given the large number of detectors 6 and the multiple rings of detectors 14, tubes-of-response 28 are typically defined (and therefore coincident events recorded) for only a portion of the pairs of detector 6.

Referring now to Figure 1 b, for example, in a PET scanner comprising rings 14 each of 256 detectors 6, 32 tubes-of response 28 are defined between each given detector and 32 opposing detectors. In the preferred embodiment, the PET scanner would therefore comprise 256 X 32, i.e. 8192 tubes-of-response 28.
Referring now to Figure 1c, in general, the detectors 6 at the ends of the tubes-of response 28 will be selected such that groups of tubes-of response run in parallel to one another, and such that tubes-of-response of different groups are not parallel. In the present example there would be 256 groups of 32 parallel tubes-of-response, with the angles between the tubes-of-response in consecutive, adjacent groups being 360/256 degrees.
Referring now to Figure 2b in addition to Figures 1 a, 1 b and 1 c, a pseudo inverse matrix P+ 30 is disclosed wherein each tube-of response 28 corresponds to a column 32 of the pseudo inverse matrix P+ 30. For an image of 64 X 64 = 4096 pixels, column 32 comprises 64 X 64 = 4096 elements 34 respectively relating a coincident event that has occurred in the associated tube-of-response to the 4096 pixels of the image. Therefore, the pseudo inverse matrix P+ 30 for a given ring 14 of detectors 6 in the preferred embodiment consists of 8192 columns and 4096 rows. The number of pixels may be safely reduced for a given image grid by not taking into account the pixels at the corners of the square image grid during the computation of the system matrix;
see 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, No. 3, pp. 1168-1175, 2000.

The index of the tube-of-response 28 of the coincident event 26 detected at time TO has a one to one correspondence with a column 32 of the pseudo inverse matrix P+ 30. The image prior to TO 36 is used as basis for the first updated image 38 which is derived by adding the image prior to TO 36 to the column 32 of the pseudo inverse matrix P+ 30 corresponding to the coincident event 26.
A second coincident event 40 is detected at time T1. The index of the tube-of response 42 of the second coincident event 40 detected at time T1 has a one to one correspondence with a second column 44 of the pseudo inverse matrix P+ 30. The updated image 38 is used as basis for the second updated image 46 which is derived by adding the updated image 38 to the second column 44 of the pseudo inverse matrix P+ 30 corresponding to the second coincident event 40.
Continuing in a similar fashion, a third coincident event 48 is detected at time T2. The index of the tube-of response 50 of this third coincident event detected at time T2 has a one to one correspondence with a third column 52 of the pseudo inverse matrix P+ 30. The second updated image 46 is used as basis for the third updated image 54 which is derived by adding the second updated image 46 to the third column 52 of the pseudo inverse matrix P+ 30 corresponding to the third coincident event 48.
Corrections of ~nro~iection data In order to provide the most accurate rendition of the tomographic image, it is appropriate to take into account physical factors and factors related to particular instruments that affect the measurement of the projection data during image reconstruction. Factors including the spatially variant system response and the model of the emission and detection processes can be taken into account when generating the system matrix. Other factors related to the specific instrument being used as well as the subject under study, such as the normalization of detector efficiency or correction for signal attenuation in tissue (e.g., photon attenuation in the case of emission tomography), can be included by multiplying the elements of regularized pseudo-inverse matrix P+ by appropriate factors, either before reconstruction or on the fly during reconstruction. The updating Equation (20) will change slightly in the case of correction on the fly:
x(S)(t)=~Cks)(t): xks)(t)=xk~t-1~+FS x per, k=1,...,M (22) where FS is the correction factor assigned to the tube-of response s.
Similarly, random coincidence events in PET could be accounted for by subtracting one column of matrix P+ from the previous image estimate:
x(S)(t)- t"kS)ltl: xks)(t)=xk~t-1~-p,~, k=1....,M (23) or, instead, by skipping the next column addition given by Equation (20) when the next coincident event is registered. In the case of correction on the fly, Equation (23) would be replaced by:
x(S)~t~=~xks)(t~: xks)(t)=xk~t-1)-FS x per, k=1,...,M (24) Hardware implementation The reconstruction of tomographic images using the above method can be implemented as a sum of two vectors representing the previous image 5 estimate and the current update from a new measured event or additional single (or individual) projection data. Additionally, it may be necessary to adjust the value of the current vector by multiplying it by a correction factor prior to addition to the previous image estimate. As discussed above, the correction factor serves to take into account any physical irregularities or irregularities in instrumentation.
Preferably, data processing is performed in an environment supporting floating point arithmetic to satisfy requirements as to precision, although in some special cases integer, or fixed point, arithmetic may be used. One major consideration which inevitably arises when attempting to carry out calculations in a fixed point environment is scaling. Scaling is necessary to maintain accuracy when a high resolution is required for signals having a wide dynamic range.
Additionally, there is the possibility of value overflow due to the limited number of bits used for the integer representation. Due to this, quantization from regions-of-interest in the reconstructed image is preferably carried out using floating point arithmetic to avoid overflow and ensure the required precision.
Real-time reconstruction subjects the processing of data to the additional constraint that the time necessary to reconstruct the image must be less than the average time between successive events. This includes the time necessary for retrieving the column of matrix P+ defined by the bin index of the current coincident event, multiplying the column by a correction factor if necessary, adding the M elements of the column to the previous image estimate and finally storing the updated image.
Depending on the manner in which projection data is being generated by the tomographic scanner, data may be acquired in one of two different ways, sequentially or randomly.
In the sequential acquisition mode, the measuring instrument sequentially supplies projection data which is measured only once in each projection bin. X-ray CT scanners and ultrasound probes are examples of tomographic scanners which fall into this category. In these cases the average time between events is determined by the scanning speed of the instrument, and is typically in the millisecond range or more. Data processing in these cases may be performed using a general-purpose computer.
In the random acquisition mode events are registered randomly in all available projection bins. SPECT and PET scanners, which operate in counting mode, are examples of tomographic scanners which fall into this category.
Assuming a peak event rate of one million counts per second, the processing of data would have to be terminated in less than one microsecond. For a relatively small 2D image (e.g., 64x64 pixels, M=4096), the requisite computational speed can merely be reached with current general-purpose computer technology by clever implementation of the algorithm. Additionally, the system matrix P+, having dimensions N x M, must also be loaded into main memory for fast access. Although extending the reconstruction to larger image sizes or 3D
imaging geometry makes such implementation impractical with the current computer technology, this should become feasible in the foreseeable future given the constant increase in both memory capacity and processing speed.
Referring now to Figure 3, a schematic diagram of a preferred image processing subsystem 56 in accordance with the present invention is disclosed.
The image processing subsystem 56 may implement a high level of parallel processing and storage capacity necessary to efficiently perform the data processing required by the real-time reconstruction of tomographic images based on the matrix pseudo-inverse method. The image processing subsystem 56 can be implemented using general purpose high-performance processors, Digital Signal Processors (DSP), programmable logic devices (e.g., Field Programmable Gate Arrays or FPGA) or application specific integrated circuits (ASIC) by persons of ordinary skill in the art.
In the preferred embodiment of Figure 3, a memory unit 58 of the image processing subsystem 56 is provided for storing the above mentioned pseudo inverse matrix P+ 30. This memory unit 58 comprises a large number of sixty-four (64) bit memory locations 60 arranged in an array of 8192 columns 62 and 4096 rows 64. As indicated in the foregoing description, each of the 8192 tubes-of response such as 28 in Figures 1 a, 1 b and 1 c of the PET scanner (not shown) is associated to a respective one of the 8192 columns of the memory unit 58.
Each column 62 comprises 4096 memory locations 60 to relate a coincident event that has occurred in the associated tube-of response to the 4096 pixels (or voxels) of the image. In operation, an index 68 associated to a tube-of-response in which a coincident event has occurred is supplied as data 82 from a PET
scanner (not shown). This index 68 is placed in an address register 84 which causes the image data from the column 62 corresponding to the tube-of-response index 68 to be placed on a data bus 66 for transfer to an optional multiplying unit 70. The image data received in the multiplying unit 70 via the data bus 66 is multiplied by a factor 72 (see correction factor Fs of Equation (23)) prior to its transfer to an adder 74. In the adder 74 the image data received from the multiplying unit 70 is used to update the image in an image buffer 76.
The contents of the image buffer 76 are made available via a display unit bus to a display unit 80.
Results A series of PET scans were performed on a high-resolution animal tomograph. The PET scanner consisted of 2 rings of 256 BGO crystals individually coupled to avalanche photodiodes. Referring now to Figure 4, the singular value spectrum of the system matrix P for the above PET scanner and an image of 64x64 pixels in 2D is shown. In this case, the matrix condition number c is equal to 4441.5.
Sequences of images demonstrating the use of real time TSVD
reconstruction are presented in both Figure 5 and Figure 6. Data was acquired using a phantom of 110 mm diameter fabricated from Lucite~. Holes having diameters of 2, 3.4, 6.7, 9.7, 13, 15.8, 20.3, 22.7 mm and located on a circumference at a distance of 28 mm from the center were machined in the phantom. Each hole was then filled with 18F-fluorodeoxyglucose, placed inside the PET scanner and the image of the acquired events reconstructed.
Referring to Figure 5, proceeding from left to right and from top to bottom each successive image is reconstructed from approximately 23800 additional counts (i.e. coincident events) in comparison to the number of counts used to reconstruct the previous image estimate. In Figure 5 the coincident events were detected randomly at all angles by a ring of detectors.
Referring now to Figure 6 an image reconstructed from the same data as that used in Figure 5 is disclosed. However, instead of proceeding sequentially in time the images are rendered using all the coincident data recorded by a detector bank, i.e. a subset of detectors, each detector in a given detector bank being selected such that pairs of detectors within the bank define parallel, or almost parallel, tubes-of-response. Proceeding from left to right and top to bottom, each successive image is reconstructed from approximately 29,750 additional counts (i.e. coincident events) in comparison to the number of counts used to reconstruct the previous image estimate, with the counts used in each successive image being detected by a eight (8) new detector banks as described above. The selection of detector banks is such that the angle of the coincident events rotated around the phantom, which allows data from a progressively increasing number of incidence angles to be included in the image reconstruction. At the outset of image reconstruction, therefore, only one or a few angles along one direction are available which gives rise to streaking artifacts appearing across the reconstructed image along the direction of incidence. Referring to the initial series of reconstructed images in Figure 6, these artifacts are quite visible. As the image reconstruction progresses with data recorded from successive detector banks being added image detail is increased. In order to achieve good results, projection angles through 180 degrees are generally required.
In both Figures 5 and 6, less than 1/3 of the singular value spectrum was used. An approach for singular value spectrum truncation based on spatial resolution analysis may be found in V. Selivanov & R. Lecomte, "Fast PET
image reconstruction based on SVD decomposition of the system matrix", ", IEEE Trans. Nucl. Sci., vol. 48, no. 3, pp. 761-767, 2001.
Discussion of application in medical imac~in4 The 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. It should be pointed out that there are other means of obtaining the inverse system matrix, for example, through direct analytical inversion of the Radon transform for individual projection data. The image of the filtered backprojected data, which represents the contribution of an individual event to the FBP reconstructed image, can be used as the corresponding column of the inverse system matrix, and the set of images of filtered backprojected data for all projections define the inverse system matrix.

Singular value spectrum truncation allows separation of the signal and noise at the reconstruction step. Index T, as discussed above, determines the trade-off between noise and resolution. Truncation of the singular spectrum is not the only way of solution regularization. Spectrum modification without 5 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 sufFciently well conditioned, no SVD pseudo-inversion is necessary and direct matrix inversion will be feasible.
TSVD reconstruction however has some drawbacks: negative values in the image estimate, streak artifacts with low-count images if the index T is lower than a certain threshold value or noise artifacts if T is higher than the threshold value. But these features (except for the noise artifacts 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 FBP
image reconstruction, the most popular method used in medical practice today.
Regardless of its drawbacks, TSVD (as well as MSVD) displays some very attractive benefits:
1 ) A 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 unnecessary since the geometry of a given system is utilized;
3) The resolution in reconstructed images may be adjusted, based on the spatial resolution analysis in reconstructed images, by varying the truncation index T (or modifying the singular value spectrum in the MSVD

case);
4) Noise amplification may be controlled by varying the truncation index T
(or modifying the singular value spectrum in the 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 the measured projection data can be reconstructed "on the fly" as soon as acquired by the scanner, storage of the measured data or intermediary calculation results, except for the image estimate being updated, is unnecessary, thus allowing for optimal 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 of the present invention. Moreover, the application of the present invention is not limited to medical imaging only, but has possible applications in other imaging techniques utilizing tomographic principles and image reconstruction from projections.

Claims (22)

1. A method for reconstructing an image, formed of an array of pixels, from a series of measured individual tomographic projection data, comprising:
creating a system matrix defining vectors relating the measured individual projection data to the pixels of the image;
associating one of the vectors of the system matrix to each measured individual projection data; and for each measured individual projection data of the series, updating the image by adding to a previous image estimate said one vector of the system matrix associated to said measured individual projection data.

2. An image recontructing method as defined in claim 1, comprising conducting said vector associating and said image updating in real-time while the individual projection data are measured.

3. An image reconstructing method as defined in claim 1, wherein creating a system matrix comprises:
factoring a system matrix defining vectors relating the measured individual projection data to the pixels of the image; and inverting the factored system matrix.

4. An image reconstructing method as defined in claim 3, wherein the factored system matrix presents the form P=UWV T , and wherein inverting the system matrix comprises pseudo-inverting the factored system matrix of the form P=UWV T to obtain a pseudo-inverse system matrix of the form P+=VW+U T, where U and V are orthogonal matrices, W is a diagonal matrix containing singular values, W+ is the inverted diagonal matrix W, and T denotes transpose of matrices V and U.

5. An image reconstructing method as defined in claim 4, further comprising truncating the pseudo-inverse system matrix to remove small singular values.

6. An image reconstructing method as defined in claim 4, further comprising modifying a singular value spectrum to reduce the effect of inversion of very small singular values.

7. An image reconstructing method as defined in claim 6, wherein modifying the singular value spectrum comprises replacing the singular values in the diagonal matrix W+ by a function of said singular values.

8. An image reconstructing method as defined in claim 1, wherein:
the measured individual projection data comprise events detected through tubes-of-response of a tomographic scanner;
the system matrix comprises a number of vectors equal to the number of tubes-of-response;
each tube-of-response is associated to a respective vector of the system matrix;
associating one of the vectors of the system matrix to each measured individual projection data comprises associating each detected event to the vector of the system matrix associated to the tube-of response through which said event has been detected; and updating the image comprises, for each detected event, adding to the previous image estimate said vector of the system matrix associated to the tube-of-response through which the event has been detected.

9. An image reconstructing method as defined in claim 1, wherein:
the measured individual projection data comprise events detected through tubes-of-response of a tomographic scanner;
the system matrix comprises a number of columns equal to the number of tubes-of-response and a number of rows equal to the number of pixels of the image;
the columns of the system matrix define respective vector-columns relating the events detected through the tubes-of response to the pixels of the image;
each tube-of-response is associated to a respective vector-column of the system matrix;
associating one of the vectors of the system matrix to each measured individual projection data comprises associating each detected event to the vector-column of the system matrix associated to the tube-of-response through which said event has been detected; and updating the image comprises, for each detected event, adding to the previous image estimate said vector-column of the system matrix associated to the tube-of-response through which the event has been detected.

10. An image reconstructing method as defined in claim 9, further comprising, prior to adding to the previous image estimate the vector-column, scaling said vector-column with a correction factor.

11. An image reconstructing method as defined in claim 1, comprising obtaining the previous image estimate by adding together all vectors of the system matrix associated to all measured individual projection data prior to measurement of the individual projection data under consideration.

12. An apparatus for reconstructing an image, formed of an array of pixels, from a series of measured individual tomographic projection data and a system matrix defining vectors relating the measured individual projection data to the pixels of the image, comprising:
a generator of a projection data index in response to each measured individual projection data;
a memory unit storing the system matrix, wherein the memory unit comprises an address input supplied with the projection data index from said generator and a data output delivering a vector of said system matrix corresponding to the projection data index;
and for updating the image for each individual projection data, an adder for adding to a previous image estimate said vector of the system matrix delivered on the data output of the memory unit.

13. An image recontructing apparatus as defined in claim 12, comprising means for operating the generator, memory unit and adder in real time while the individual projection data are measured.

14. An image reconstructing apparatus as defined in claim 12, wherein the system matrix is a pseudo-inverse system matrix of the form P+=VW+U T, which results from a peuso inversion of a system matrix having the form P=UWV T, where U and V are orthogonal matrices and W is a diagonal matrix containing singular values, W+ is the inverted diagonal matrix W, and T
denotes transpose of matrices V and U.

15. An image reconstructing apparatus as defined in claim 14, further comprising means for truncating the pseudo-inverse system matrix to remove small singular values.

16. An image reconstructing apparatus as defined in claim 14, further comprising means for modifying a singular value spectrum to reduce the effect of inversion of very small singular values.

17. An image reconstructing apparatus as defined in claim 16, wherein said means for modifying the singular value spectrum comprises means for replacing the singular values of the diagonal matrix W+ by a function of said singular values.

18. An image reconstructing apparatus as defined in claim 12, wherein:
the measured individual projection data comprise events detected through tubes-of response of a tomographic scanner;
the memory unit comprises a number of columns of memory locations equal to the number of tubes-of-response and a number of rows of memory locations equal to the number of pixels of the image;
the columns of the memory unit store respective vector-columns relating the events detected through the tubes-of-response to the pixels of the image;
each tube-of-response is associated to a respective column of the memory unit; and the adder has a first input receiving the previous image estimate and a second input receiving the vector-column stored in the column of the memory unit associated to the tube-of-response through which the event has been detected.

19. An image reconstructing apparatus as defined in claim 18, further comprising a multiplier which, prior to adding to the previous image estimate the vector-column, multiplies said vector-column with a correction factor.

20. An image reconstructing apparatus as defined in claim 12, comprising means for obtaining the previous image estimate by adding together all vectors of the system matrix associated to all measured individual projection data prior to measurement of the individual projection data under consideration.

21. An apparatus for reconstructing an image, formed of an array of pixels, from a series of measured individual tomographic projection data and a system matrix defining vectors relating the measured individual projection data to the pixels of the image, comprising:
means for associating one of the vectors of the system matrix to each measured individual projection data; and means for updating the image for each measured individual projection data of the series, said updating means comprising an adder for adding to a previous image estimate said one vector of the isystem matrix associated to said measured individual projection data.

22. An image reconstructing apparatus as recited in claim 21, comprising means for operating the associating means and the updating means in real time while the individual projection data of the series are measured.