Abstract:
The invention concerns a method which consists in processing in a central unit ( 4 ) the digitized image ( 1 ) obtained by sensing noisy radiation. The image ( 1 ) is considered as a table of pixel intensity values, which is broken down into p elementary tables of n pixels so as to subsequently order them into a processing table with p lines of n columns; applying to said table a factorial analysis method, to deduce therefrom the n significant factors; then reconstituting a reconstructed processing table taking into account the most significant factors, and in deducing therefrom a reconstituted image wherein the high-frequency noise is reduced, maintaining a satisfactory contrast.

Description:
TECHNICAL FIELD OF THE INVENTION  
       [0001]     The present invention relates to methods and devices for filtering digitized images, in particular to methods and devices for eliminating noise in a digitized image.  
         [0002]     Some images are very noisy, which reduces their legibility. This is the case, for example, of images obtained from low-level signals, in which the wanted signals are severely disturbed by interferences and other physical phenomena that degrade signal transmission.  
         [0003]     This is also the case of certain images in the field of medical imaging, such as scintigraphy images.  
         [0004]     In nuclear medicine, the scintigraphy technique consists in injecting patients with biological molecules marked by radioactive isotopes that emit gamma photons. The gamma photons are detected by a gamma camera.  
         [0005]     The emitter particles injected into the human body emit gamma photons at random. The signals picked up depend, on the one hand, on the quantity of emissive particles in a given region and, on the other hand, on the emission conditions for each particle. In a scintigraphy image, which is an image of the distribution of gamma photons from an organism, each picture element or pixel is an integer and is the number of particles, i.e. the number of photons detected facing it. It is found that the distribution of these integer values around a mean value is a Poissonian distribution.  
         [0006]     The Poissonian distribution is expressed by the formula:  
         Pr   ⁡     (     X   =   k     )       =         μ   k     ⁢           ⁢     e     -   μ           k   !           
 
         [0007]     Pr (X=k) is the probability of the pixel X taking the value k. μ is the mean value of the distribution.  
         [0008]      FIG. 1  shows one example of a Poissonian distribution in which the mean μ=5.  
         [0009]     Another instance of a Poissonian distribution is the distribution of X-ray photons such as those used in tomodensitometry (i.e. in an X-ray scanner), which visualizes the anatomical structures of an organism using differences in the attenuation of X-rays by biological tissues.  
         [0010]     In all cases, the images are affected by Poissonian noise, resulting from the random character of the emissions, and the relative magnitude of this noise is inversely proportional to the number of events detected. This is because the relative error of the measurement is given by the formula:
 
Relative error=[(variance) ( X )] 1/2 /mean( X )={square root}μ/μ=1/{square root}μ
 
         [0011]     The relative error therefore varies in the opposite sense to the number of events detected.  
         [0012]     To improve the legibility of an image, it is therefore only natural to reduce the noise, and to this end to increase the number of events detected.  
         [0013]     A first way to do this is to acquire images over a longer time period. However, this may be difficult if the subject is moving, for example in radiology because of respiratory and cardiac movements. Moreover, increasing the image acquisition time ties up the sensing apparatus for longer, which raises serious problems in scheduling medical imaging.  
         [0014]     A second way to increase the number of photons detected is to increase the flux of particles. In medical imaging, this implies increasing the radioactive dose injected or the flux of X-ray photons. However, this commensurately increases the dose of radiation to which the patient is exposed, which contradicts the principles of radio protection.  
         [0015]     Finally, using more-efficient sensing devices may be envisaged, but this immediately gives rise to cost problems.  
         [0016]     To illustrate the effect on image quality of the increased number of photons detected,  FIG. 2  reproduces a real scintigraphy image of a subject, produced from signals received over three different periods: the left-hand image corresponds to reception for one minute, the center image corresponds to reception for five minutes, and the right-hand image corresponds to reception for ten minutes. The quality of the right-hand image is very much higher.  
         [0017]     Poissonian noise clearly consists in a random variation between one pixel and the adjacent pixels, variation that is not related to the real difference in the concentration of emissive particles between the region facing the pixel concerned and the regions facing the adjacent pixels. This noise is therefore high-frequency noise, i.e. noise introducing sudden variations from one pixel to another in the sequence of the pixels of an image.  
         [0018]     As it is often very difficult to increase the number of photons detected, the invention proposes to reduce the high-frequency noise contained in an image by using a particular filtering device.  
         [0019]     Filtering images, with the aim of reducing the effect of noise, has already been proposed. A very large number of different filters has been proposed.  
         [0020]     The simplest and best known type of filter replaces the value of each pixel by the mean of the values of the adjacent pixels. A Gaussian filter has also been proposed, which imparts a Gaussian spatial distribution weight to the adjacent pixels, and replaces the central pixel with a linear combination of the adjacent pixels.  
         [0021]     The document WO 01/72032 describes one example, which defines particular coefficients for applying filtering and, where possible, at the same time preserves a certain contrast in the image.  
         [0022]     The document also states the problem of all filters, whether weighted or not: all filters improve the signal to noise ratio, but at the cost of severely degrading the spatial resolution and the contrast of the image. This is because a filter calculating a weighted mean of a plurality of adjacent pixels obviously reduces sudden variations of intensity between adjacent pixels, and therefore reduces the effects of contrast.  
         [0023]     By way of example,  FIG. 3  depicts the processing of a real scintigraphy image by two filters: the left-hand view is an unfiltered raw image; the center view is the same image filtered by a 3×3 weighted filter; the right-hand view is the same image filtered by a 3×3 median filter, which replaces each pixel with the median of the adjacent pixels. It can be seen that the contours, in the filtered images, are not very sharp, as a result of the reduction in spatial resolution and contrast.  
         [0024]     Filtering images by multivariate statistical analysis methods applied to a table of numbers each expressing the degree of luminosity of a corresponding pixel in a digitized image has also been proposed.  
         [0025]     Thus, the document XP 002212651 AEBERSOLD, STADELMANN, ROUVIERE describes processing electron microscope images. This processing is applied by factorial analysis of correspondences to a table consisting of the pixels of the digitized image. The analysis is applied to the entire digitized image, seeking in the entire image the correspondence factors most representative of the image, i.e. those that correspond to the highest eigenvalues, and then reconstructing the image on the basis of only these representative factors. Although in theory the method achieves high quality filtering of a given image, in practice it is not possible to adapt the filtering to the structure of the image to be filtered, as a result of which results are disappointing.  
         [0026]     The document XP 001074312 HANNEQUIN, LIEHN, VALEYRE proposes applying factorial analysis of correspondences to a series of entire scintigraphy images, and suggests using the likelihood ratio test to determine the correspondence factors to be used to reconstruct the series of images. The method is again applied to the totality of the images. The results of such filtering are slightly improved, but remain disappointing in respect of certain image structures to be filtered.  
         [0027]     The same problems are encountered in the documents XP 008007773 and XP 008007666, which apply statistical processing to the totality of an image.  
         [0028]     There is therefore a need for a more efficient filtering device, that eliminates high-frequency noise without significantly reducing the resolution or the contrast of the image.  
         [0029]     One problem common to all the filtering devices is the difficulty of adapting the filtering process to the structure of the image to be filtered. This is because the images do not all have the same structure, the structure differing as a function of the subject observed. It appears that although filters may generally be well adapted to a particular image, they are ill adapted to very different images. Tests have been carried out to adapt the filters as a function of the variations of the pixels, but these tests are not satisfactory.  
       SUMMARY OF THE INVENTION  
       [0030]     Thus an object of the present invention is to enable efficient filtering, to produce images of high quality, using means that are particularly inexpensive.  
         [0031]     In one particularly advantageous embodiment, the invention therefore proposes to adapt the filtering process not only to the general structure of an image, but also, and above all else, to the individual structure of basic portions of the image to be processed. Clearly this may significantly improve the contrast of the image.  
         [0032]     Furthermore, according to the invention, this adaptation is preferably automatic, thereby minimizing the loss of information resulting from the filtering.  
         [0033]     To achieve the above and other objects, the invention provides a method of processing a digitized image consisting of a table T of numbers x (ij)  each expressing the degree of luminosity of a corresponding pixel of rank i, j (for example, the number of particles detected in the corresponding pixel), the method comprising the reduction of high-frequency noise (for example Poissonian statistical noise) by the following steps:  
         [0034]     a) decomposing the table T into a continuous series of p basic tables of the same dimensions each having n pixels,  
         [0035]     b) ordering the data from the series of basic tables into a processing table X of p rows and n columns, each row i being formed of the ordered sequence of the pixels from the i th  basic table,  
         [0036]     d) effecting a multivariate statistical analysis on the processing table X, considering the n columns as variables, to extract therefrom n representative factors,  
         [0037]     f) generating a reconstituted processing table XR of numbers xr (i, j)  using only the first q representative factors and restoring the absolute degree of luminosity, then generating a reconstituted table TR constituting the reconstituted digitized image in which the high-frequency noise has been reduced in this way,  
         [0038]     during the step f), the reconstituted processing table XR is reconstituted by reconstructing each row i independently of the others, taking into account only the factors having a meaningful weight with the row i.  
         [0039]     In the present description and in the claims, the expression “multivariate statistical analysis processing” refers generally to all statistical analysis methods such as in particular analysis into main components with or without normalization, factor analysis with or without normalization, and factor analysis of correspondences.  
         [0040]     Thanks to row-to-row reconstruction, the invention adapts the filtering process to each subset of the image, which considerably improves the result obtained.  
         [0041]     One advantageous embodiment of the method uses a method of factorial analysis of correspondences, whereby:  
         [0042]     d1) the transposed matrix XN T  is calculated,  
         [0043]     d2) the square matrix that is the product of the normalized matrix XN and the transposed matrix XN T  is calculated,  
         [0044]     d3) the square matrix XN T XN is diagonalized to extract therefrom n eigenvectors u k  (with k from 1 to n) associated with n eigenvalues vp k ,  
         [0045]     d4) the coordinates of the p rows of the normalized matrix XN are calculated on the n eigenvectors,  
         [0046]     d6) the coordinates of the n columns are calculated on the n eigenvectors.  
         [0047]     Advantageously: 
        during a step d5) the squared cosines of the rows on the n factor axes are calculated,     the squared cosines are used to test the weight of the representative factors.        
 
         [0050]     In this case, after the step b), a normalization step c) can advantageously be provided for obtaining a normalized matrix XN in which each element xn ij  of row i and column j is weighted by a transform TN using the mean of the values of the elements of the row i and the mean of the elements of the column j, during which normalization step:  
         [0051]     c1) the sum f i  of each row of the processing table X is calculated,  
         [0052]     c2) the sum f j  of each column of the processing table X is calculated,  
         [0053]     c3) the total sum f tot  of the table is calculated,  
         [0054]     c4) a normalization transform is used to replace each element x ij  of the table by the normalized value equal to x ij  divided by the product of the square roots of f i  and f j .  
         [0055]     In one advantageous embodiment, the invention proposes an auto-adaptive reconstruction method, to eliminate as much high-frequency noise as possible without affecting the quality of the image. The means employed to this end consist in retaining only the variance of the signal and eliminating the variance of the noise. To this end, according to the invention, the reconstruction of a row of the reconstituted table XR is stopped as soon as a sufficient number of factors have been used for the variance of the reconstructed row to be greater than the variance of the original signal from which the estimated variance of the noise has been subtracted.  
         [0056]     The test applies particularly well to processing an image affected by Poissonian noise, for example. In this case, one property of Poisson&#39;s law is that its mean is equal to its variance:
 
variance ( X )=mean ( X )=μ
 
         [0057]     As a result of this a good estimate of the variance of the noise is obtained by calculating the mean of the signals from row i.  
         [0058]     In another aspect, the invention proposes a device for processing digitized images, comprising a memory, a calculation unit, an input-output device for receiving data constituting a digitized image to be processed, display and/or printing means for viewing the processed image, and a program stored in memory and adapted to execute the method as defined hereinabove.  
         [0059]     The invention applies in particular to a medical imaging installation comprising a device of the above kind. 
     
    
     BRIEF DESCRIPTION OF THE DRAWINGS  
       [0060]     Other objects, features and advantages of the present invention will emerge from the following description of particular embodiments of the invention given with reference to the appended drawings, in which:  
         [0061]      FIG. 1  depicts a Poissonian distribution with a mean equal to 5;  
         [0062]      FIG. 2  depicts the improvement of an image by increasing the duration of observation of a Poissonian phenomenon;  
         [0063]      FIG. 3  depicts results obtained by prior art weighted or median filters;  
         [0064]      FIG. 4  depicts the result that may be obtained with one embodiment of a filter according to the invention;  
         [0065]      FIG. 5  depicts the main steps of a filtering method according to the present invention;  
         [0066]      FIG. 6  depicts an offsetting principle used in the method of the invention to limit edge effects; and  
         [0067]      FIG. 7  depicts diagrammatically an imaging installation incorporating an image processing device according to the invention. 
     
    
     DESCRIPTION OF THE PREFERRED EMBODIMENT  
       [0068]     Consider first  FIG. 7 , which depicts diagrammatically a medical imaging installation comprising a gamma ray sensor  1  such as a gamma camera moving in two directions in front of the subject  2  to be observed, and capturing gamma rays  3  emitted by radio-emissive particles injected beforehand into the body of the patient  2 . The gamma ray sensor  1  sends to a calculation unit  4  the sequence of photon image signals received on each basic region or pixel of the gamma ray sensor  1 , the calculation unit  4  storing in a memory  5  the number of photons for each pixel, corresponding to the intensity of the pixel. The memory  5  therefore contains a digitized image consisting of a table T of numbers x (i, j) each expressing the number of photons detected (or degree of luminosity) of a pixel from the row i and the column j of the observed region  1 . According to the invention, the installation further comprises a program stored in the memory  5  and driving the calculation unit  4  to filter the image digitized in the above manner and to produce on a display or printing device  6  a filtered image of high quality from which high-frequency noise has been extracted.  
         [0069]     An image filtering method according to one embodiment of the present invention is described next, with reference to  FIG. 5 .  
         [0070]     Note the table T already referred to and constituting the digitized image. In practice, digitized images contain a large number of pixels. To simplify the explanation, a square image is considered that comprises 8×8 pixels, each pixel being depicted by a small square.  
         [0071]     The first operation a) of the noise reduction method according to the invention consists in decomposing the table T into a continuous series of p basic tables of the same size, each comprising n pixels. In the example depicted in  FIG. 5 , four basic tables T 1 , T 2 , T 3  and T 4 , each having 16 pixels, are considered.  
         [0072]     Then, in a step b), the data from the sequence of basic tables T 1 -T 4  is arranged in a processing table X of p rows and n columns, each row i being formed of the ordered sequence of the pixels of the basic table of rank i. Thus pixels 1 to 16 from the table T 1 , stored in order, are repeated in the first row of the table X. Similarly, the pixels stored in order from the table T 2  are repeated in the second row of the table X, and so on. Thus, in this example, the processing table X has four rows each of 16 columns.  
         [0073]     In the  FIG. 5  example, the table T is decomposed into four square basic tables T 1 -T 4 . However, without departing from the scope of the invention, the table T may be divided into a sequence of rectangular tables of the same size.  
         [0074]     Then, in a step c), the processing table X is normalized to obtain a normalized matrix Xn in which each element xn ij  in row i and column j is weighted by a transform using the mean of the values of the elements of the row i and the mean of the elements of the column j. One example of a normalization transform TN is described later.  
         [0075]     Then, in the step d), AFC statistical factor analysis processing is applied to the normalized matrix XN, taking the j columns as the variables, to extract therefrom n representative orthogonal factors.  
         [0076]     Factor analysis of correspondences is used, for example, in which:  
         [0077]     d1) the transposed matrix XN T  is calculated,  
         [0078]     d2) the square matrix that is the product of the normalized matrix XN and the transposed matrix XN T  is calculated,  
         [0079]     d3) the square matrix XN T XN is diagonalized to extract therefrom n eigenvectors u k  (with k from 1 to n) associated with n eigenvalues vp k ,  
         [0080]     d4) the coordinates of the p rows of the normalized matrix XN are calculated on the n eigenvectors,  
         [0081]     d5) the squared cosines of the p rows of the n factor axes are calculated,  
         [0082]     d6) the coordinates of the n columns are calculated on the n eigenvectors.  
         [0083]     During a step e), the n factors are then stored in decreasing order as a function of their respective weights.  
         [0084]     During a step f), a reconstituted processing table XR of numbers xr (i, j) is generated using only the first q representative factors and restoring the original values of the pixels by means of a transform that is the inverse of the normalization transform TN.  
         [0085]     Finally, the reconstituted table TR is generated, constituting the reconstituted digitized image, in which high-frequency noise, for example Poissonian statistical noise, has been reduced.  
         [0086]     In one embodiment, the step c) may apply normalization by means of the following steps:  
         [0087]     c1) the sum f i  of each row of the processing table X is calculated,  
         [0088]     c2) the sum f j  of each column of the processing table X is calculated,  
         [0089]     c3) the total sum f tot  of the table is calculated,  
         [0090]     c4) a normalization transform is used to replace each element x ij  of the table by the normalized value equal to x ij  divided by the product of the square roots of f i  and f j .  
         [0091]     In practice, the step d4) of calculating the coordinates of the p rows of the normalized matrix XN on the n eigenvectors may be executed by calculating coordinate c k  (i) of the row i on the axis k generated by the eigenvector U k  using the formula:  
           c   k     ⁡     (   i   )       =       f   tot     1   /   2       ⁢       ∑     j   =   1     n     ⁢         x   ij         (       f   i     ⁢           ⁢     f   j       )       1   /   2         ⁢           ⁢       u   k     ⁡     (   j   )                 
 
         [0092]     in which u k (j) is the j th  coordinate of the eigenvector u k .  
         [0093]     The squared cosine may be calculated, according to the step d5), using the formula:
 
cos 2   k ( i )=c k ( i ) c k ( i )
 
         [0094]     The coordinates of the n columns of the n eigenvectors may be calculated during a step d6) from the coordinate d k (j) of the column j on the eigenvector of axis k using the formula:  
           d   k     ⁡     (   j   )       =       1     vp   k     1   /   2         ⁢       ∑     i   =   1     p     ⁢         x   ij       f   j       ⁢           ⁢       c   k     ⁡     (   i   )                 
 
         [0095]     in which vp k  is the eigenvalue associated with the eigenvector u k .  
         [0096]     During the step e), the n factors can advantageously be classed as a function of their squared cosine by applying the above formulas.  
         [0097]     According to the invention, the reconstitution of the reconstituted processing table xr is effected independently row by row, taking into account only the q factors having the maximum squared cosine for the row i.  
         [0098]     Assuming that the first q factors are taken into consideration, the reconstructed value xr ij  (q) of the element of the reconstituted processing table XR in row i and column j is calculated from the formula:  
           xr   ij     ⁡     (   q   )       =           f   i     ⁢           ⁢     f   j         f   tot       ⁢       ∑     k   =   1     q     ⁢           c   k     ⁡     (   i   )       ⁢           ⁢       d   k     ⁡     (   j   )           vp   k     1   /   2                 
 
         [0099]     From the reconstituted processing table XR, in  FIG. 5 , the reconstituted table TR is reconstructed row by row, the first row of the reconstituted processing table XR constituting the pixels of the first basic table TR 1 , and so on.  
         [0100]     According to the invention, the aim is to automate the adaptation of the filtering device to the content of the image. This automation is effected for each region of the image corresponding to one of the basic tables T 1  to T 4 . To this end the reconstituted table XR is reconstructed row by row. Reconstructed values xr ij  of a row i of elements of the reconstituted table XR are calculated step-by-step: 
        the values of the elements xr ij  are calculated successively for increasing values q of the number of factors taken into account,     the residual variance of the row i is calculated each time,     the residual variance is compared to the estimated variance of the noise to be reduced,        
 
         [0104]     and the calculation for the row i is stopped when the residual variance of the row i is no longer statistically greater than the estimated variance of the noise of the row i in the starting image, thus obtaining an estimated final image Im_final.  
         [0105]     In practice, the residual variance var_res(q) is calculated as the difference between the initial variance of the row i of the processing table X and the reconstituted variance var_rec(q) or variance of the row i of the reconstituted processing table XR:
 
var_res( q )=var_ini−var_rec( q )
 
         [0106]     The test of comparing the residual variance and the estimated variance of the noise may advantageously be effected by:  
         [0107]     a) calculating the variable t using the formula:
 
 t =(Var_noise) xhi ( ddl )/ ddl 
 
         [0108]     in which xhi (ddl) is the value given by the χ 2  table for a risk of 5% and a number ddl of degrees of freedom,  
         [0109]     ddl is the number of degrees of freedom, with
 
 ddl=n−q −1
 
         [0110]     q is the number of factors taken into account,  
         [0111]     b) stopping reconstruction when the residual variable Var_res(q) is less than t.  
         [0112]     If the method is applied to processing an image subject to Poissonian noise, the estimated variance of the noise of the row i is taken as equal to the average of the elements x ij  of the row i of the processing table X.  
         [0113]     To reduce the effect of noise on the results, the procedure described above is repeated several times on the same image, shifting the decomposition into basic tables by one pixel each time.  FIG. 5  depicts the first procedure for an offset of 0 in x and 0 in y.  FIG. 6  depicts the second procedure for an offset of 1 pixel in x and 0 pixels in y: the basic table T′ 1  is offset by 1 pixel toward the right in the table T. For a decomposition into 4×4 squares, for example, the procedure is executed 16 times, with offsets in x from 0 to 3 and offsets in y from 0 to 3. The final image (Im_final), estimated without noise, is the mean of the 16 resulting images.  
         [0114]     This mean can take into account the number of times each pixel of the image is actually included in the AFC factor analysis of correspondences processing, in order not to cause edge effects. Another advantage of repetition is circumventing geometrical artifacts that may appear because of the decomposition into basic rectangles.  
         [0115]     An advantageous embodiment of the invention further chooses the level of elimination of noise, subtracting from the original image only a portion of the noise image, to obtain a so-called reduced image im_reduced. The image of the noise im_noise is the original image subtracted from the final image. Accordingly, a reduced image Im_reduced is produced in which the noise is partially eliminated by the following steps: 
        calculating the image of the noise Im_noise as the difference between the original image and the final image:
 
Im_noise=Im_original−Im_final
    partially adding the image of the noise to the final image:
 
Im_reduced=Im_final+βIm_noise
       
 
         [0118]     α may range from 0 (total elimination of the noise) to 1 (no elimination of the noise).  
         [0119]      FIG. 4  depicts the result of filtering in accordance with the invention: the left-hand view is the original image; the middle view is the final image with the optimum filtering; the right-hand view is the reduced image, with α=⅓.  
         [0120]     The present invention is not limited to the embodiments explicitly described and encompasses variants and generalizations thereof within the scope of the following claims.