Image processing systems and methods
FIELD OF THE INVENTION 
     The present invention relates generally to image processing systems and methods for producing and modifying images used for example in clinical studies of the heart, brain, lungs, kidneys and other organs, and more particularly to systems and methods for optimally enhancing said images. 
     In particular, the inventive system improves the signal-to-noise ratio in images which are statistically limited, and results in images comparable to those acquired by prior art methods using substantially longer time periods or substantially more radiation. 
     BACKGROUND OF THE INVENTION 
     A number of medical diagnostic imaging systems known in the art are examples of systems which would be improved by implementing the basic concept of the present invention. 
     In general, the prior art systems utilize detectors to sense radiation such as X or γ flux interacting with the detectors and to convert each quantum of X- or γ-rays to individual electrical output pulses or to convert integrated inputs to proportional electrical outputs. Location-computing circuitry calculates the X and Y coordinates of each interaction and thus identifies the elemental area of the detector or array of detectors within which an interaction (event) occurs to a given accuracy. The systems further process and display images produced from the detected data. For example, the present invention may be implemented as a data processing means for a Gamma Camera (Scintillation Camera) of the Anger type. Such cameras are widely used in Nuclear Medicine as radiation detectors for determining the distribution of a radioactive substance in a body. The camera converts images formed by gamma radiation into a related visible image. See, for example, U.S. Pat. No. 3,011,057 in which typical Gamma Camera apparatus is disclosed for collecting information from patients after the administration of radio-nuclides through inhalation, ingestion or injection. 
     The present invention can also be used in Computerized Axial Tomography (CAT) and ultrasound imaging, as well as in non-diagnostic systems such as starlight intensifiers, radar systems and so on. Apparatus for performing Computerized Axial Tomography, for example, produces a representation of the absorption coefficients of penetrating radiation from a source of X-rays at a plurality of elemental locations distributed over a cross-sectional slice of a body under investigation. Such apparatus well-known to the persons skilled in the art is disclosed, for instance, in the instance, in the U.S. Pat. No. 3,778,614. 
     Regardless of the nature of the body being examined, and the type of radiation, only those beams incident on a radiation detector in accordance with the present invention contribute to the outputs of the detector. 
     As used herein the &#34;actual physical data&#34; that is being measured, processed and used, means either a distribution (spatial or projected) of absorption density, such as tissue density in the case of X-rays, a distribution of radioisotope concentration in the case of γ-rays, a distribution of tissue interfaces in the case of ultrasound detection, or a distribution of radar reflectors, optical sources or reflectors and so on, depending on the case. 
     The actual physical data obtained from Gamma Camera or CAT systems includes significant amounts of &#34;statistical noise&#34;. 
     The &#34;statistical noise&#34; is the random variation between the &#34;true&#34; value, which would have been achieved by imaging the object with infinite radiation and the actual measured values which also depend on the statistical variability of finite radioactive or absorption events. Conventional data processing techniques, for example those producing and displaying images obtained from Gamma Camera or CAT systems, include arrangements for reducing the statistical noise. More particularly, smoothing of different types, filters of different types and correctors of different types are provided, all of which are based on some &#34;global&#34; assumption such as the known aperture of the system or the known resolution limit below which everything is considered as &#34;statistical noise&#34;. These &#34;global&#34; assumptions, as used in the prior art, have often tended to degrade the images by either actually increasing the noise or reducing the resolution, because of variations between the local case and the global assumptions. Thus, there is a long felt need for means other than prior art arrangements for reducing the statistical noise. 
     BRIEF DESCRIPTION OF THE INVENTION 
     It is an object of the present invention to provide image-optimizing techniques generally and particularly for use with medical diagnostic imaging radiation detectors. The image-optimizing techniques, in accordance with one aspect of the present invention, are incorporated in digital image-processing systems used with medical radiation apparatus. They are used to calculate estimators of the actual physical data being measured in a way that maximizes the likelihood that the estimators are correct. 
     In this way the image optimizing techniques help to reduce the amount of statistical &#34;noise&#34; in the image. 
     The present technique, although incorporating whatever knowledge we have of the aperture or resolution, does not treat these known factors as global constants; instead, the corrections to the actual physical data are performed during the suggested optimizing process, using regional and semilocal data. 
     The method used to reduce the amount of statistical &#34;noise&#34; in the image by the present invention is based in part on an &#34;estimation&#34; technique developed by Stein and others for sampling and polling. (See &#34;Stein&#39;s Paradox in Statistics&#34; by Bradley Efron and Carl Morris, Scientific American, May 1977, pp 119-127.) Stein&#39;s theory, as incorporated into applied statistics, concerns the use of the averages of observed quantities to estimate unobservable quantities. In accordance with James and Stein an estimator for unobservable quantities may be determined by the relation Z=Y+C (Y-Y), where Y is the average of a single set of data, Y is the grand average of averages and C is a constant defined as a &#34;shrinking factor&#34;. 
     Essential differences between the Stein-James estimator and the image optimizing method and system used in accordance with the present invention as well as other features and advantages of the invention will become more apparent as the invention becomes better understood from the following detailed description made in conjunction with the accompanying drawings, wherein: 
    
    
     BRIEF DESCRIPTION OF THE DRAWINGS 
     FIG. 1 is a block diagram representation of a medical radiation detector system in accordance with the present invention; 
     FIG. 2 is a plan schematic view of a detector of the system of FIG. 1 in accordance with the present invention; 
     FIG. 3 is a graphical representation of the accuracy of locating an event; 
     FIG. 4 is a flow diagram of the present invention; and 
     FIG. 5 is a schematic representation of the definition of regions around different points in the matrix. 
    
    
     GENERAL DESCRIPTION 
     FIG. 1 at 10 shows a medical diagnostic imaging system for carrying out the image optimizing method in accordance with the present invention. The head or front end is designated by numeral 21. It generally comprises a source of X-ray or γ-flux, and detector means 11 such as at least one scintillation crystal. Transducer means such as an array of photomultipliers 14 is shown mounted proximate to the scintillator. The photomultipliers convert the light output of the scintillators to electrical signals. The electrical signal outputs of the photomultipliers are amplified in matrix amplifiers 16. 
     Position computing circuitry 17, for computing x and y coordinates of an event and thus identifying the elemental area within which an event occurs, are provided. The circuitry 17 receives location information and the energy signals for the particular locations from the matrix amplifier 16. The computed, x,y location of each event is at the output of circuitry 17. 
     Means are provided for processing the signals received from circuitry 17. More particularly processing and display equipment 26 is provided. The equipment includes digital computing circuitry 22 for counting the electrical pulses connected with the arrival of X-ray or -photons at the said detector means and for storing the original event counts (N) in a two-dimensional matrix memory 23. A digital processor 24 transforms the number (N) of events at given x and y locations to processed data (F); for example, by multiplying the said number (N) by a given arbitrary normalization factor G. The obtained processed data (F) is then stored in the matrix at the location corresponding to the x,y location of the event. 
     The matrix memory 23 in the present invention may reside either in the internal memory of the processor or on some external memory like a magnetic tape, a disk or a floppy disk. 
     Although the description is drawn to an embodiment in which the matrix (23) has two dimensions, i and j, representing the x and y coordinates, it is evident to those skilled in the art that the same methods and techniques are also valid for matrices having one dimension i or three dimensions (i,j,k) where the third dimension represents such characteristic as either the third space dimension or the time dimension. Also, although the description is drawn to x- or γ-ray detection, the methods may also apply to other energy means used for producing images such as (but not limited to) phonons, infrared radiation, radio wave radiation and microwave radiation. 
     An arithmetic unit 25 shown in FIG. 1 may be incorporated in the present invention into the digital processor 24, used for modifying the stored data (F) and replacing the stored data by the modified data, providing a less noisy image on display unit 27. The less noisy image is provided in the present invention by following the image-optimizing techniques fully illustrated in FIG. 4. 
     This image-optimizing phase accepts as inputs a normalization factor (G), two weights W1 and W2 associated with the stored data (F), and the matrix of the stored data (Fi) corresponding to element i defined by position coordinates xi and yi. 
     The normalization factor (G) can be either the known normalization factor or an arbitrary constant or a factor calculated from the stored data (F). The stored data (F) in the conventional radiation-measuring process is generally retained in a matrix of columns and rows. The matrix locations are specifically related to elemental areas at corresponding locations on the detector. The data stored in the matrix is generally accepted as being an estimator of some property of the organ or tissue being examined. For example, in a Gamma Camera the stored number of counts is taken to represent the number of radioisotope emissions striking the detector along the line of sight. However, the generally accepted procedure for working with known Gamma Cameras and X-ray tomographs produces images whose quality is limited by the total event counts (statistics limited). Due to the randomly irregular nature of radioactive emission and radiation absorption the same amount of radioisotope emission or radiated organ absorption observed for the same length of time will result in different event counts each time. It is usually accepted that the average number from many such observations is the &#34;true&#34; estimator, and the differences between different measurements and this value are considered statistical noise. 
     Therefore a present need exists for reducing such noise by using radiation-measuring equipment offering techniques which provide a new estimator value closer to the &#34;true&#34; estimator. 
     The reduction of such noise is presently best achieved by using more photons; such as by either increasing the radiation dose or the examination duration. Both of these prior art solutions have severe drawbacks. Other methods provide a new estimator obtained from the acquired data without increasing the radiation dose or the examination time. The new estimator value is obtained by an examination of the data in neighbouring line-of-sight locations and reduces the number of photons required per line-of-sight. The data in the neighbouring line-of-sight locations or an estimator based on such data can be used in place of measured values, if the temporal &#34;noise&#34; is approximately equal in magnitude to the spatial &#34;noise&#34;. This assumption is justified if due notice is taken of the true morphological variation in the amount of radiosotope emission from point to point. 
     An essential part of the invention is the fitting of a two-dimensional surface representative of this variation in the measured quantity of events from point to point to the stored measured data (F) around the point of interest. The residual variation is representative of the statistical noise (both temporal and spatial) and can be used as a measure of the success of the fit. For regions that are small enough, two dimensional 2nd order surfaces usually provide a close fit of the morphology. For larger regions, higher order functions are necessary. The fitted surfaces are relatively insensitive to the local &#34;noise&#34; and become less sensitive with larger regions. By using optimizing techniques for bivariate probabilities, it can be concluded that the best estimator that can be derived from the data (F) stored in the matrix and the fitted surface is: C*=αC+(1-α)C; where: 
     C*--is the new calculated estimator; 
     C--is a regional estimator as derived from the fit of the two-dimensional surface to the measured data; 
     C--is a local estimator from the stored data, usually taken as the stored data itself at the point of interest; C=F(X,Y)=Fi; 
     α--a localized weight function depending on C, C and V--the regional variability of the processed data (F) around the said two-dimensional surface. 
     The arithmetic unit 5 accomplishes the optimizing process, accepting the input mentioned above. From this input the arithmetic unit 5 proceeds by using the processed data (F) stored in a given matrix element, together with the processed data (F) stored in neighbouring matrix elements, for instance, in accordance with the following steps: 
     1. Choosing a region around each point corresponding to a matrix element, the region size being determined by the assumed accuracy of the identification of the X and Y coordinates, said point serves as a center point for the region; 
     2. Fitting a two-dimensional surface to the processed data (F) of all the points (i.e. matrix elements) inside said region; 
     3. Calculating the regional estimator C at the said point, for example by taking C to be the value of the said surface at the center point of the region; 
     4. Obtaining the regional spread V, for example, by taking the variability of the processed data around the said surface in the said region; 
     5. Computing the regional weight α and the local weight (1-α), wherein α=(W 2  C)/(W 1  V+W 2  C); where W1 and W2 depend on assumptions made on the nature of the processed data (F), and for the case of event counting are calculated by the formulae: 
     
         W.sub.1 =(2C/G)/(C/G+C); 
    
     
         W.sub.2 =2C/(C/G+C) 
    
      and G is the said normalization factor; 
     6. Obtaining the new estimator C* by multiplying the regional estimator (C) by the regional weight, the local estimator (C) of the given point, which can be Fi by the local weight (1-α), and summing the two products. The above weights are calculated to maximize the likelihood of the new estimator being correct: C*=αC+(1-α)C. 
     At the end, means are provided for displaying the improved image obtained by the above process and means are provided for replacing the stored data (F) by the improved results or for storing both. Thus, for obtaining this new &#34;true&#34;estimator, the data which was spread around each point was used to correct the value (data) which was measured at the point. The techniques and methods, in accordance with the present invention, optimize the probability of correctness and thus simulate the effects of obtaining a larger event count, which would necessitate using undesirable methods such as lengthening the acquisition time or enlarging the dose of the radio-active pharmaceutical or the X-ray flux. 
     In addition it should be noted that the image-optimizing techniques for obtaining the &#34;true&#34; estimator substantially differ from Stein&#39;s method and other methods derived from Stein&#39;s methods. These differences arise from the following factors: Stein&#39;s method and the methods following it were either one-dimensional or had no well-ordered two-dimensional arrangement, such as demographical data. 
     Digital images, on the other hand, have a well-ordered two-dimensional arrangement and the correlation between close neighbours falls off with distance in a fairly location-independent way. This fact makes it possible to apply a region of constant shape and size around each point, within which we assume all points to be highly correlated with the center point (see FIG. 5). 
     The nature of the physical data is such that all values of the event-counts, processed data, fitted data, optimized data and so on, are non-negative over the whole image. This constraint causes the local weight (1-α) and the regional weight (α) to be in the region of between 0 and 1. Such constraints were not an inherent part of Stein&#39;s method. 
     The regional fit of the surface to the data utilizes only the data that could have come from the same aperture as the data at the point being processed. Demographical data have no such physical aperture, therefore a fit is usually done with a global variable or at least a regional variable of a region which is very large and having a shape which depends on the case. From this follows that whereas the arithmetic unit (5) performs a fit over a different, although similarly shaped region for each point being processed (see FIG. 5). Stein&#39;s method after defining the region uses the same regional estimator for all points within the region. 
     The same sort of consideration forced Stein and his followers to fit the data to a constant which is usually the regional average. The present invention, on the other hand, takes into account the morphology of the object under examination by fitting the data inside each region to a more general two-dimensional surface. 
     The nature of the physical data used in imaging is such that it is statistically limited. The physical data thus comes with a criterion for its reliability, namely, the mean relative error of the data at each point is equal to the inverse of the square root of the original number of counts on which this data is based. This fact is incorporated in the weights used in the arithmetic unit (25) resulting in regional and local weights much different from those used by Stein which have no such inherent reliability measure. The nature of the fit also indicates the use of the relative variability of the data around the fitted surface as the measure of reliability of the regional estimator. This is consistent with the reliability measure used for the local estimator and replaces the standard deviation used by Stein and his followers. 
     As a result of all the above points the regional weight (α) obtained herein is independent of the local estimator (C) and also of the local &#34;noise&#34; (√C), unlike other weighting methods where α depends on both local and regional estimators and on the reliability measures of both of them. For the case of event counts, the above regional weight (α) is given by: 
     
         α=C/(V/G+C) 
    
     The above methods of the image optimizer of the present invention elevate the performance of the medical diagnostic imaging systems. 
     While the invention has been particularly shown and described with reference to several embodiments from nuclear medicine, it will be evident to those skilled in the art that it can be implemented in industrial, military, astronomical, scientific, and other imaging systems including other medical imaging systems, without departing from the spirit and the scope of the invention. It being understood that all these uses are encompassed within the scope of the invention as determined by the accompanying claims.