Abstract:
The invention relates to the field of digital image processing and can find use in suppression of noise in digital images, formed by high-energy radiation, including X-ray radiation. Specifically, the invention relates to a method for suppression of noise in digital x-ray images. The objective of the invention is to provide a method for noise suppression in digital X-ray images that has extended functionality, specifically, a method that makes it possible to reduce residual noise level and amount of artifacts in the form of discontinuities, directed along local orientation of object borders in textured image areas, to reduce residual LF noise level, and to eliminate over-smoothing (excessive filtering) of fine details. The technical innovation achieved is the improvement of digital image processing quality.

Description:
RELATED APPLICATIONS 
       [0001]    This application claims priority to Eurasian Patent Application No. EA 201200255, filed Feb. 22, 2012, which is incorporated herein by reference in its entirety. 
       FIELD OF THE INVENTION 
       [0002]    The invention relates to the field of digital image processing and can find use in solving problems of noise suppression in digital images, formed by high-energy radiation, including X-ray radiation. Specifically, the invention relates to a method for suppression of noise in digital x-ray images. 
       BACKGROUND OF THE INVENTION 
       [0003]    In current clinical practice, various image processing methods are used to correct and facilitate interpretation of digital X-ray images, including: sharpening, enhancing of anatomical structure imaging, automatic targeting of range of interest etc. Some of the most important methods for X-ray image enhancement include noise suppression methods. Noise suppression ensures better visual perception and distinctiveness of low-contrast details and improves quality of execution of other image processing algorithms. Therefore, one of high priority tasks in modern digital radiology is to develop noise filtering methods, capable of preserving diagnostically relevant information and providing high-speed processing at significant levels of noise, dependent on signal intensity. 
         [0004]    One can identify two major problems, faced by a researcher who undertakes to develop a digital X-ray image filtering algorithm. The first problem refers to estimation of image noise level, which essentially depends on signal intensity. The second problem refers to development of a filtering method that would ensure the right balance between filtering quality and filtering speed. A closer look should be taken at each of these problems. 
         [0005]    The problem of estimation of noise, dependent on signal intensity, can be summarized as follows. Practically any qualitative method for digital image noise filtering involves use of data on amplitude and spectral distribution of noise. In fairly wide range of noise suppression methods, noise variance value is utilized as a parameter. Publications by various authors often use a model, which addresses to digital image noise as a stationary random process with normal amplitude distribution and spectral distribution of white noise. Also, studies of non-stationary noise suppression problem are underway, to be more specific, with regards to noise with signal-dependent variance, and it has been shown that the more accurate is the known noise variance value, the more effective is the noise suppression process [Bosco et. al., Noise Reduction for Cfa Image Sensors Exploiting Hvs Behavior, Sensors 9(3), 1692-1713, 2009; Foi, Practical Denoising of Clipped or Overexposed Noisy Images, Proc. 16th European Signal Process. Conf., EUSIPCO, 2008]. 
         [0006]    The nature of digital X-ray image noise should be examined. A detector measures intensity of attenuated (by passing through the object examined) X-ray radiation. During imaging each detector cell accumulates the average of  N  electrons due to absorption of photons. The quantity of accumulated electrons N can be simulated as a Poisson distributed random variable: 
         [0000]    
       
         
           
             
               
                 P 
                  
                 
                   ( 
                   
                     N 
                     = 
                     n 
                   
                   ) 
                 
               
               = 
               
                 
                   exp 
                    
                   
                     ( 
                     
                       - 
                       
                         N 
                         _ 
                       
                     
                     ) 
                   
                 
                  
                 
                     
                 
                  
                 
                   
                     
                       N 
                       _ 
                     
                     n 
                   
                   
                     n 
                     ! 
                   
                 
               
             
             , 
             
               n 
               ≥ 
               0. 
             
           
         
       
     
         [0007]    Accidental fluctuations of the number of absorbed photons are called photon noise. In modern detectors, photon noise is the primary source of noise. Additional sources include detector system noises, such as: reading noise, thermal noise, amplifier noise, quantization distortion etc. Cumulative effect of such noises is simulated as a Gaussian distributed random variable [Bosco et. al.,  Noise Reduction for Cfa Image Sensors Exploiting Hvs Behavior, Sensors  9(3), 1692-1713, 2009; Foi et. al.,  Practical poissonian - gaussian noise modeling and fitting for single - image raw - data, Image Processing, IEEE Transactions , TIP-03364-2007-FINAL]. According to a widely used model, in linear electronic circuits variance of photon noise and additional noises is linearly dependent on useful signal value [Yane,  Digital image processing , Moscow,  Technosphera,  2007] 
         [0000]      σ 2 ( u ( p ))= au ( p )+ b,   (1)
 
         [0008]    where u(p) represents signal intensity level in pixel p=(i, j) of an image, and σ 2 (u(p)) represents noise variance dependence on signal intensity. 
         [0009]    Therefore digital image noise is linearly dependent on signal intensity. There is a considerable number of publications on the problem of estimation of signal-dependent noise [Bosco et. al., Noise Reduction for Cfa Image Sensors Exploiting Hvs Behavior, Sensors 9(3), p. 1692-1713, 2009; 
         [0010]    Foi et. al., Practical poissonian-gaussian noise modeling and fitting for single-image raw-data, Image Processing, IEEE TRANSACTIONS, TIP-03364-2007-FINAL. 
         [0011]    Förstner, Image Preprocessing for Feature Extraction in Digital Intensity, Color and Range images, In Springer Lecture Notes on Earth Sciences, 1998; 
         [0012]    Hensel et. al., Robust and Fast Estimation of Signal-Dependent Noise in Medical X-Ray Image Sequences, Springer, Proceedings BVM 2006, Hamburg, Germany, March 19-21, p. 46-50, 2006; 
         [0013]    Liu et. al., Automatic estimation and removal of noise from a single image, IEEE Transactions on Pattern Analysis and Machine Intelligence, p. 1-35, October, 2006; 
         [0014]    Salmeri et. al.,  Signal - dependent Noise Characterization for Mammographic Images Denoising.  16 th  IMEKO TC4 Symposium (IMEKOTC4 &#39;08), Florence, Italy, 2008; 
         [0015]    Waegli, Investigations into the Noise Characteristics of Digitized Aerial Images, In: Int. Arch. For Photogr. And Remote Sensing, Vol. 32-2, 1998]. For example, the following study [Hensel et. al., Robust and Fast Estimation of Signal-Dependent Noise in Medical X-Ray Image Sequences, Springer, 2006] refers to establishing a distribution-free method for frame-by-frame estimation of noise in series of digital X-ray images, specifically emphasizing development of an algorithm, suitable for real-time execution. The study [Foi et. al., Practical poissonian-gaussian noise modeling and fitting for single-image raw-data, Image Processing, IEEE Transactions on 17, 1737-1754, 2008] describes a two-parameter approach to digital image noise estimation. Noise of original digital image (as captured by detector and not yet subjected to non-linear transformations such as gamma correction etc.) is simulated based on a random variable, which represents an additive to the signal: 
         [0000]        u   n ( p )= u ( p )+σ 2 ( u ( p )) n ( p ),  (2)
 
         [0016]    where u n (p) represents signal intensity level in pixel p of the noisy image observed, and σ 2 (u(p)) represents noise variance dependence on intensity of a signal of type (1), and n(p)εN(0,1) represents a standard normal random variable. The above study advances a method for tracing noise variance simulation curves, which takes into account sensor nonlinearities, causing underexposure and overexposure, i.e. delinearity at the edges of dynamic range. 
         [0017]    The second problem refers to necessity to achieve balance between filtering quality and high-speed filtering requirement. In modern X-ray units, the necessity for rapid high-resolution image processing is somewhat contradictory to commitment to high filtering quality. By now, there has been developed and put to use an extensive mathematical apparatus, related to digital image noise suppression challenge. In general, virtually any noise suppression method can be assigned to one of two classes of filtering methods: 
         [0018]    filtering in image transformation space; 
         [0019]    filtering in original image space. 
         [0020]    The first group of methods (filtering in transformation space) refers to conversion of an original image by applying some transformation. Such transformation is followed by processing transformation factors based on a certain rule, e.g. by applying an appropriate nonlinear operator. This kind of processing is aimed at reduction of noise energy and, possibly, enhancing image sharpness. The final stage of the filtering process is image reconstruction by inverse transformation. This group includes a wide range of methods based on: Fourier transform, wavelet transform [Strang et. al.,  Wavelets and Filter Banks , Wellesley-Cambridge Press.; Vetterli et. al.,  Wavelets and Subband Coding , Prentice Hall, Englewood Cliffs, N.J., 1995], complex wavelet transform [Kingsbury,  The dual - tree complex wavelet transform: A new technique for shift invariance and directional filters, In Proc.  8 th IEEE DSP Workshop , Utah, Aug. 9-12, 1998, paper no. 86; Selesnick,  The double - density dual - tree discrete wavelet transform. IEEE Trans. Signal Processing , vol. 52, no. 5, pp. 1304-1314, May 2004], curvelet transform [Candes et. al.,  Curvelets, multiresolution representation, and scaling laws, In SPIE conference on Signal and Image Processing: Wavelet Applications in Signal and Image Processing VIII , San Diego, 2000] etc. Execution speed and filtering quality in this class of methods are determined by transformation type and factor processing method used. 
         [0021]    The second group of methods (filtering in original image space) includes equally extensive list of approaches, including: anisotropic diffusion [Weickert,  Anisotropic Diffusion in Image Processing . Tuebner Stuttgart, 1998], variation based algorithms [Rudin et. al.,  Nonlinear total variation based noise removal algorithms, Physica D,  60:259-268, 1992], bilateral filtering [Tomasi et. al.,  Bilateral Filtering for Gray and Color Images, In Proc.  6 th Int. Conf Computer Vision , New Delhi, India, 1998, pp. 839-846], non-local means filtering [Buades et. al.,  On image denoising methods, SIAM Multiscale Modeling and Simulation,  4(2):490-530, 2005], regression filtering [Takeda et. al.,  Kernel regression for image processing and reconstruction, IEEE Transactions on Image Processing , vol. 16, no. 2, pp. 349-366, February 2007], optimal spatial adaptation filtering [Kervrann et. al.,  Optimal spatial adaptation for patch - based image denoising. IEEE Transactions on Image Processing , vol. 15, no. 10, October 2006] etc. As a rule, such methods reduce to searching for similarity in pixels or groups of pixels (blocks or patches) within the image itself, possibly based on geometric data (such as local orientation of image structures), and application of some pattern for modification (equalization) of brightness of such pixels, whereby modification rate is dependent on degree of similarity of pixels (pixel blocks). 
       SUMMARY OF THE INVENTION 
       [0022]    The present invention proposes a method for noise filtering in digital X-ray images that includes two stages. The first stage involves development of estimation of signal-dependent noise. The second stage involves noise suppression and includes a number of additional substeps, such as: combined bilateral non-local filtering, pyramid reduction of residual low-frequency noise and over-smoothing correction. 
         [0023]    The closest to the digital X-ray image noise estimation method disclosed herein is the method described in the following publication [Hensel et. al., Robust and Fast Estimation of Signal-Dependent Noise in Medical X-Ray Image Sequences, Springer, 2006], according to which the following operations are sufficient to build an estimation of signal-dependent noise based on an original image: 
         [0024]    Estimate useful signal by means of low-frequency filtering of the original image and calculate difference between the original image and its estimation, thus forming a noise image; 
         [0025]    Use some appropriate method to remove pixel values of the noise image that are associated with abrupt changes of the original image (borders and single hot pixels); 
         [0026]    Divide intensity range of the estimation image into intervals and, with regard to each interval, accumulate values of noise image pixels, associated with noise estimation image pixels; 
         [0027]    Calculate noise variance in each intensity interval based on values of noise image pixels, accumulated in such interval. 
         [0028]    The present invention utilizes principles of noise estimation, outlined in the above publication. The most significant distinctions include: (a) deletion of values of noise image pixels, associated with abrupt changes of original image, by means of morphological enhancement of the values of noise image pixels, associated with original image borders; and (b) use of robust local linear approximation of noise variance interval estimations to create a tabulated function that describes noise-signal dependency; and (c) based on estimation image and the obtained tabulated function, calculation of noise map in the form of an image that represents a pixel-by-pixel estimation of noise variance in the original digital image. 
         [0029]    The nearest equivalent of the method for noise suppression in digital X-ray images, used in the present invention, is non-local means filtering [Buades et. al., On image denoising methods, SIAM Multiscale Modeling and Simulation, 4(2):490-530, 2005]. Non-local means filter shows good results when used to suppress noise in digital images. This filter uses area averaging to estimate current pixel signal intensity: 
         [0000]    
       
         
           
             
               
                 
                   
                     
                       u 
                       ^ 
                     
                      
                     
                       ( 
                       p 
                       ) 
                     
                   
                   = 
                   
                     
                       ∑ 
                       
                         q 
                         ∈ 
                         
                           S 
                            
                           
                             ( 
                             
                               p 
                               , 
                               t 
                             
                             ) 
                           
                         
                       
                     
                      
                     
                       
                         w 
                          
                         
                           ( 
                           
                             p 
                             , 
                             q 
                           
                           ) 
                         
                       
                        
                       
                         
                           u 
                            
                           
                             ( 
                             q 
                             ) 
                           
                         
                         / 
                         
                           
                             ∑ 
                             
                               q 
                               ∈ 
                               
                                 S 
                                  
                                 
                                   ( 
                                   
                                     p 
                                     , 
                                     t 
                                   
                                   ) 
                                 
                               
                             
                           
                            
                           
                             
                               w 
                                
                               
                                 ( 
                                 
                                   p 
                                   , 
                                   q 
                                 
                                 ) 
                               
                             
                             . 
                           
                         
                       
                     
                   
                 
               
               
                 
                   ( 
                   3 
                   ) 
                 
               
             
           
         
       
     
         [0030]    where S(p,t) represents an area of radius t, centered at the current pixel p (i.e. area of 2·t+1×2·t+1 pixels in size); w(p,q)=F(∥U(p, f)−U(q, f)∥ 2,a ,σ,k) represents weight of averaging, and F represents a decreasing function; U(p, f), U(q, f) represent pixel blocks of radius f, centered at pixels p and q, respectively; 
         [0000]    
       
         
           
             
               
                  
                 
                   
                     U 
                      
                     
                       ( 
                       
                         p 
                         , 
                         f 
                       
                       ) 
                     
                   
                   - 
                   
                     U 
                      
                     
                       ( 
                       
                         q 
                         , 
                         f 
                       
                       ) 
                     
                   
                 
                  
               
               
                 2 
                 , 
                 a 
               
             
             = 
             
               
                 
                   ∑ 
                   
                     h 
                     ∈ 
                     
                       S 
                        
                       
                         ( 
                         
                           0 
                           , 
                           f 
                         
                         ) 
                       
                     
                   
                 
                  
                 
                   
                     G 
                      
                     
                       ( 
                       
                         h 
                         , 
                         a 
                       
                       ) 
                     
                   
                   · 
                   
                     
                       ( 
                       
                         
                           u 
                            
                           
                             ( 
                             
                               p 
                               + 
                               h 
                             
                             ) 
                           
                         
                         - 
                         
                           u 
                            
                           
                             ( 
                             
                               q 
                               + 
                               h 
                             
                             ) 
                           
                         
                       
                       ) 
                     
                     2 
                   
                 
               
             
           
         
       
     
         [0000]    represents Euclidean distance between pixel blocks, weighted by Gaussian function G(h,a)=exp(−h·h T /a 2 ); σ represents RMS error of image noise; k represents a parameter that defines averaging rate; u(q) represents pixel intensity in the original image, and û(p) represents resultant estimation of intensity of the current pixel. 
         [0031]    In standard interpretation, e.g. as described in [Buades et. al., Image Denoising by Non-Local Averaging, Proc. of IEEE Int. Conf. on Acoustics, Speech, and Signal Processing, vol. 2, pp. 25-28, March 2005], weight function is defined as 
         [0000]        F (∥ U ( p,f )− U ( q,f )∥ 2,a   ,σ,k )=exp(−∥ U ( p,f )− U ( q,f )∥ 2,a   2   /k   2 σ 2 );
 
         [0032]    and noise is regarded as an additive of homoscedasticity, i.e. noisy image model is defined by type (2) formula, where σ 2 (u(p))=σ 2 =const for any pixel p. Basic non-local filtering sequence, outlined above, had served as a basis for a wide range of studies, related to various aspects of improvement, acceleration and practical application of non-local means filter [Bilcu et. al.,  Fast non - local means for image de - noising, In Digital Photography III Proc. of SPIE IS &amp; T , Russel A. Martin, Jeffrey M. DiCarlo, and Nitin Sampat, Eds., February 2007, vol. 6502; Dauwe et. al.,  A Fast Non - Local Image Denoising Algorithm, In Proc. SPIE Electronic Imaging , vol. 6812, San Jose, USA, January 2008; Mahmoudi et. al.,  Fast image and video denoising via nonlocal means of similar neighborhoods, IEEE Signal Processing Letters,  12(12).839-842, 2005; Peter et. al.,  Robust estimation approach for nonlocal - means denoising based on structurally similar patches, Int. J. Open Problems Compt. Math ., vol. 2, no. 2, June 2009; Wang et. al.,  Fast non - local algorithm for image denoising, In Proc. of IEEE International Conference on Image Processing  ( ICIP ), 2006, pp. 1429-1432]. 
         [0033]    Non-local means filter is notable for high computational complexity. One of possible ways to speed up operation of non-local means filter was described in [Buades et. al.,  A review of image denoising algorithms, with a new one, Multiscale Modeling and Simulation  ( SIAM interdisciplinary journal ), Vol. 4 (2), pp. 490-530, 2005], where it is proposed to use vectorial non-local mean. 
         [0000]    In that algorithm version, image I(p) is partitioned into the M of blocks U(p m , f), m=1,M, 
         [0000]    
       
         
           
             
               
                 I 
                  
                 
                   ( 
                   p 
                   ) 
                 
               
               = 
               
                 
                   
                     ⋂ 
                     M 
                   
                   
                     m 
                     = 
                     1 
                   
                 
                  
                 
                   U 
                    
                   
                     ( 
                     
                       
                         p 
                         m 
                       
                       , 
                       f 
                     
                     ) 
                   
                 
               
             
             , 
           
         
       
     
         [0000]    where p m  represents a central pixel of partition block. Partition blocks are assumed to overlap. Estimation of intensity levels of pixels of the current partition block is made based on the following formula 
         [0000]    
       
         
           
             
               
                 
                   
                     
                       U 
                       ^ 
                     
                      
                     
                       ( 
                       
                         
                           p 
                           m 
                         
                         , 
                         f 
                       
                       ) 
                     
                   
                   = 
                   
                     
                       ∑ 
                       
                         q 
                         ∈ 
                         
                           S 
                            
                           
                             ( 
                             
                               
                                 p 
                                 m 
                               
                               , 
                               t 
                             
                             ) 
                           
                         
                       
                     
                      
                     
                       
                         w 
                          
                         
                           ( 
                           
                             
                               p 
                               m 
                             
                             , 
                             q 
                           
                           ) 
                         
                       
                        
                       
                         
                           U 
                            
                           
                             ( 
                             
                               q 
                               , 
                               f 
                             
                             ) 
                           
                         
                         / 
                         
                           
                             ∑ 
                             
                               q 
                               ∈ 
                               
                                 S 
                                  
                                 
                                   ( 
                                   
                                     
                                       p 
                                       m 
                                     
                                     , 
                                     t 
                                   
                                   ) 
                                 
                               
                             
                           
                            
                           
                             
                               w 
                                
                               
                                 ( 
                                 
                                   
                                     p 
                                     m 
                                   
                                   , 
                                   q 
                                 
                                 ) 
                               
                             
                             . 
                           
                         
                       
                     
                   
                 
               
               
                 
                   ( 
                   4 
                   ) 
                 
               
             
           
         
       
     
         [0034]    where w(p m ,q)=exp(−∥U(p m , f)−U(q, f)∥ 2   2 /k 2 σ 2 ) represents blocks averaging weight; 
         [0000]    
       
         
           
             
               
                  
                 
                   
                     U 
                      
                     
                       ( 
                       
                         
                           p 
                           m 
                         
                         , 
                         f 
                       
                       ) 
                     
                   
                   - 
                   
                     U 
                      
                     
                       ( 
                       
                         q 
                         , 
                         f 
                       
                       ) 
                     
                   
                 
                  
               
               2 
               2 
             
             = 
             
               
                 ∑ 
                 
                   h 
                   ∈ 
                   
                     S 
                      
                     
                       ( 
                       
                         0 
                         , 
                         f 
                       
                       ) 
                     
                   
                 
               
                
               
                 
                   ( 
                   
                     
                       u 
                        
                       
                         ( 
                         
                           
                             p 
                             m 
                           
                           + 
                           h 
                         
                         ) 
                       
                     
                     - 
                     
                       u 
                        
                       
                         ( 
                         
                           q 
                           + 
                           h 
                         
                         ) 
                       
                     
                   
                   ) 
                 
                 2 
               
             
           
         
       
     
         [0000]    represents squared Euclidean distance between current pixel block and comparison pixel block; and Û(p m ,f) represents the resultant estimation of intensity levels of pixels in current partition block. Final block estimations are integrated with averaged overlapping of blocks, thus forming filtered output image. Using vector-based process allows significant filter speed up, however, with a certain loss of processed image quality. 
         [0035]    The outlined basic filtering process (3),(4) has a number of downsides: 
         [0036]    High level of residual noise and specific artifacts in the form of discontinuities, directed along local orientation of object borders in textured image areas; 
         [0037]    Presence of residual LF noise; 
         [0038]    Excessive filtering (over-smoothing) of small details. 
         [0039]    The objective of the present invention is to provide a method for noise suppression in digital X-ray images that has extended functionality, specifically, a method that makes it possible to reduce residual noise level and amount of artifacts in the form of discontinuities, directed along local orientation of object borders in textured image areas, to reduce residual LF noise level, and to eliminate over-smoothing (excessive filtering) of fine details. 
         [0040]    In terms of the claimed method for digital X-ray image noise suppression, comprising: 
         [0041]    original image acquisition; 
         [0042]    estimation of signal-dependent noise variance; 
         [0043]    combined bilateral non-local filtering; 
         [0044]    suppression of residual low-frequency noise; and 
         [0045]    over-smoothing correction, 
         [0046]    the objective, i.e. improvement of digital X-ray image processing quality, is accomplished by: 
         [0047]    morphological elimination of values of noise image pixels, associated with borders in the original image, at signal-dependent noise variance estimation stage; 
         [0048]    using robust piecewise linear approximation of interval estimations of noise variance to obtain a tabulated function of noise-signal dependency; 
         [0049]    calculating a noise map in the form of pixel-by-pixel estimation of noise variance in the original image; 
         [0050]    averaging by combining non-local means filter and bilateral filter weights at combined bilateral non-local filtering stage; 
         [0051]    using a pyramid method with residual low-frequency noise filtering by a bilateral filter at residual low-frequency noise reduction stage; and 
         [0052]    obtaining a result image and correction of over-smoothing by combining the final image and the original image based on their structural and pixel similarity. 
         [0053]    In an optional embodiment of the present invention, at combined bilateral non-local filtering stage, averaging weights are calculated based on an image, smoothed by applying linear LF filter. 
         [0054]    Another optional embodiment of the present invention involves application of combined weight function in the form of fixed weight fractions of non-local means filter and bilateral filter at combined non-local bilateral filtering stage; at the final stage outputs of these filters are combined. 
         [0055]    Another optional embodiment of the present invention involves application of separable version of bilateral filter at combined bilateral non-local filtering stage, involving separation of non-local means filtering and bilateral filtering. 
     
    
     
       BRIEF DESCRIPTION OF THE DRAWINGS 
         [0056]    The claimed technical solution, as well as potential for its practical implementation and accomplishment of its practical objective, outlined above, are illustrated in  FIGS. 1 through 5 . 
           [0057]      FIG. 1  shows X-ray imaging device 
           [0058]      FIG. 2  shows X-ray image of a head 
           [0059]      FIG. 3  shows plotted tabulated function  5  (solid line) that defines relationship between standard noise deviation and signal intensity with reference to the image in  FIG. 2   
           [0060]      FIG. 4  shows a map of standard noise deviation of the image in  FIG. 2   
           [0061]      FIG. 5  A shows effectiveness of the proposed method 
           [0062]      FIG. 5  B shows effectiveness of the proposed method 
           [0063]      FIG. 5  C shows effectiveness of the proposed method 
           [0064]      FIG. 5  D shows effectiveness of the proposed method 
           [0065]      FIG. 5  E shows effectiveness of the proposed method 
       
    
    
     DETAILED DESCRIPTION OF THE PREFERRED EMBODIMENTS 
       [0066]    X-ray images can be produced, for example, by a device shown in  FIG. 1 . It includes X-ray tube  1 , emitting X-ray beam  2 . X-ray beam  2  enters detector  3 . Detector  3  includes scintillation screen (not shown) and matrix array camera (not shown). Scintillation screen is optically connected to the active area surface of the matrix array camera. X-ray beam  2  hits detector  3 , and scintillation screen converts it into visible light, which is further converted by detector sensors into a digital image. 
         [0067]    A filter, proposed by the present invention, is used to perform 4-stage noise suppression. Each of these stages is examined in detail below. 
         [0068]    Stage 1. 
         [0069]    Estimation of noise variance, dependent on signal intensity. At this stage, original image is first estimated by applying linear low-frequency image filtering [Hensel et. al., Robust and Fast Estimation of Signal-Dependent Noise in Medical X-Ray Image Sequences, Springer, 2006]. To meet strict filtering speed requirements, it is reasonable to use simple linear filtering (e.g. by applying binomial filter). The resulting smoothed image is used to obtain a noise image, which essentially represents the difference between the original and the filtered images. 
         [0070]    Image estimation by applying simple image filters is not perfect, as image edges tend to become over-smoothed. As a result, using difference between the original image and the smoothed image leads to forming a noise image, which includes, in addition to noise pixels in smoothed areas, a certain amount of pixels, associated with abrupt changes (i.e. borders of anatomical structures; such pixels are hereafter referred to as non-noise pixels). Such pixels can considerably distort the estimated variance value and should be excluded from calculations. This can be accomplished by various methods [Foi et. al., Practical poissonian-gaussian noise modeling and fitting for single-image raw-data. Image Processing, IEEE Transactions on 17, 1737-1754 (October 2008), Salmeri et. al., Signal-dependent Noise Characterization for Mammographic Images Denoising. IMEKO TC4 Symposium (IMEKOTC4 &#39;08), Firenze, Italy, September 2008], which essentially boil down to thresholding of smoothed derivatives of the original image, where threshold value is determined by local estimation of signal/noise relationship. In image areas with large number of details such estimation usually proves to be inadequate. Therefore, as far as rejection of borders is concerned, this invention proposes a simpler approach to border enhancement that does not require calculation of derived and local estimates of standard deviation. Such morphological approach to non-noise pixels rejection can be summarized as follows: 
         [0071]    The noise image is divided into two components—binary images of positive and negative changes; 
         [0072]    To enhance areas, associated with borders, resulting images are subjected to morphological operations of erosion, followed by dilatation; 
         [0073]    Processed images are combined to receive a single binary image, i.e. a map of original image borders. 
         [0074]    To better preserve fine structures, morphological erosion and dilatation are performed by small size masks (2×2 window). 
         [0075]    Calculation of interval estimations of noise variance involves identifying minimum and maximum intensity of estimation image (intensity range edges), followed by selection of partitioning interval, e.g., equal to 32 gray gradations. After that, for each image pixel it is necessary to identify an interval that pixel&#39;s value pertains to, and to use the relevant noise image pixel value to calculate noise variance estimation within that interval (this involves eliminating pixels that are associated with borders). In calculating interval estimation of noise variance, various formulas may be applied, such as standard unbiased estimator or robust estimation, based on median absolute deviation formula [Hensel et. al., Robust and Fast Estimation of Signal-Dependent Noise in Medical X-Ray Image Sequences, Springer, 2006; Foi et. al., Practical poissonian-gaussian noise modeling and fitting for single-image raw-data, Image Processing, IEEE Transactions on 17, 1737-1754, 2008]. This procedure results in a tabulated function that defines noise variance dependence on signal intensity. 
         [0076]    Inaccuracies in generation of borders map may cause serious errors in calculation of variance interval estimations. Therefore such estimations are refined by applying iterative outlier removal technique [Hensel et. al., Robust and Fast Estimation of Signal-Dependent Noise in Medical X-Ray Image Sequences, Springer, 2006]: specifically, this involves iterative elimination in each interval of noise pixels values that exceed in absolute magnitude a threshold, equal to three standard noise deviations, followed by recalculation of noise variance estimation in that interval. 
         [0077]    Once calculation of noise variance interval estimations is completed, estimation of noise variance dependence on signal intensity is performed. Using a parametric estimator, it is basically possible to develop a noise model parameters estimation (1), using some of available methods (e.g. least-squares method, likelyhood function minimization or directed optimization). It is also possible to account for sensor nonlinearities [Foi et. al.,  Practical poissonian - gaussian noise modeling and fitting for single - image raw - data, Image Processing, IEEE Transactions on  17, 1737-1754, 2008]. However, as was noted in [Hensel et. al.,  Robust and Fast Estimation of Signal - Dependent Noise in Medical X - Ray Image Sequences , Springer, 2006], for a number of reasons, development of parametric model, able to adequately define noise variance dependence on signal intensity, can pose serious difficulties. Such reasons might include sensor operational nonlinearities, or nonlinear original image preprocessing techniques (such as logarithmation). Therefore, from the point of view of implementation convenience and application versatility, the optimum approach is the one that involves development of non-parametric noise-signal dependency estimation based on estimations of interval variance. 
         [0078]    The present invention uses non-parametric approach to determination of required relationship, which involves development of interpolating tabulated function based on identified interval estimations of noise variance. Such tabulated function is formed based on robust local linear approximation of interval estimations of noise variance. Application of robust methods allows additional reduction of outlier effects (gross errors in interval variance estimations), while locality of approximation ensures iteration of complex run of a curve, describing dependence of noise on signal intensity. Therefore the resultant tabulated function of each original image intensity associates with noise variance estimation. The function of table entry points can be performed by, e.g., intensities of the estimation image. For a diagram of such single tabulated function (solid line)  5 , plotted for the image shown in  FIG. 2 , alongside with interval estimations of RMS noise deviation  6  (circles), see  FIG. 3 . 
         [0079]    In practice, parametric noise estimation may include an approach that involves a conversion process, stabilizing noise variance in the original image [Starck et. al., Image Processing and Data Analysis:  The Multiscale Approach . Cambridge University Press, 1998, Yane,  Digital image processing , Moscow,  Technosphera,  2007]. In this case, the problem of filtering noise, dependent on signal intensity, is reduced to the challenge of suppression of additive signal-independent noise with preset variance properties. The present invention uses non-parametric noise estimation, therefore it is suggested that the following approach is used in noise map development: building, based on estimation image and interpolating table, a noise map, i.e. an image where each pixel estimates root-mean-square deviation of noise in associated pixel of the original image. The noise map provides pixel-by-pixel noise estimation at precision level, sufficient for practical applications. For a noise map, developed based on image in  FIG. 2 , see  FIG. 4 . This example illustrates performance of noise estimation algorithm in overexposure (i.e. sensor saturation/clutter) conditions that can be easily observed around the head (light halo, pos. 4 ) on  FIG. 2 . 
         [0080]    Stage 2. Combined Filtering. 
         [0081]    When standard filtering process (3) is used, as the block size increases, filtering quality, measured in peak values of signal/noise relationship, improves. In addition to quality improvement, in case of application of vectorized non-local mean (4), filter execution speed increases along with pixel block enlargement. However, in highly textured image areas and along abrupt intensity variations (borders), block enlargement leads to relatively high level of residual noise and specific artifacts in the form of discontinuities, directed along local orientation of object borders. Presence of residual noise and artifacts in textured areas can be attributed to decrease of noise suppression rate, which is caused by reduction of the quantity of such blocks. High level of residual noise can also be observed in smooth areas of processed image. The present invention relates to development of combined noise filtering method, based on combination of bilateral filter and non-local means filter, that should improve downsides of the standard method. Such method is hereafter referred to as combined filtering algorithm, or simply combined filter. Combined bilateral non-local filtering is based on the following principle. 
         [0082]    Combined Bilateral Non-Local Filtering Principle. 
         [0083]    In the course of averaging, weight of each pixel is calculated based on its proximity to the current pixel, using two metrics: bilateral (brightness) and non-local (structural). If the current pixel and the comparison pixel are close from the point of view of non-local means filter, i.e. their blocks are structurally similar, domination of structural weight should be observed. Otherwise, if blocks are not similar, yet current pixel and comparison pixel brightness levels are close, this leads to enhancement of the effect of brightness weight. To reduce residual noise, averaging weights should be calculated based on smoothed (estimation) image, which can be produced by applying a simple linear low-frequency filter to the original image. 
         [0084]    Mathematically, this principle leads to the following modification of standard filtering process (3): 
         [0000]    
       
         
           
             
               
                 
                   
                     
                       u 
                       ^ 
                     
                      
                     
                       ( 
                       p 
                       ) 
                     
                   
                   = 
                   
                     
                       ∑ 
                       
                         q 
                         ∈ 
                         
                           S 
                            
                           
                             ( 
                             
                               p 
                               , 
                               t 
                             
                             ) 
                           
                         
                       
                     
                      
                     
                       
                         w 
                          
                         
                           ( 
                           
                             p 
                             , 
                             q 
                           
                           ) 
                         
                       
                        
                       
                         
                           u 
                            
                           
                             ( 
                             q 
                             ) 
                           
                         
                         / 
                         
                           
                             ∑ 
                             
                               q 
                               ∈ 
                               
                                 S 
                                  
                                 
                                   ( 
                                   
                                     p 
                                     , 
                                     t 
                                   
                                   ) 
                                 
                               
                             
                           
                            
                           
                             
                               w 
                                
                               
                                 ( 
                                 
                                   p 
                                   , 
                                   q 
                                 
                                 ) 
                               
                             
                             . 
                           
                         
                       
                     
                   
                 
               
               
                 
                   ( 
                   5 
                   ) 
                 
               
             
           
         
       
     
         [0085]    Here, w(p,q)=V(b(p,q),n(p,q)) represents a function that defines averaging weight of combined filter, selected in such a way as to satisfy the combined filtering principle; b(p,q)=G(q−p,a)·F((ū(p)−ū(q)), σ(ū(p)),k b ,h b ) represents the weight of bilateral filter, and G(s,a)=exp(−s·s T /a 2 ) represents Gaussian function; n(p,q)=F(∥Ū(p,f)−Ū(q,f)∥ 2 ,σ(ū(p)),k n ,h n ) represents the weight of non-local means filter; ū(p), Ū(p, f) represents the respective values of pixel and pixel block in estimation image; and σ(ū(p)) represents the value of RMS error of estimation image noise, dependent on signal intensity. To determine similarity of pixels brightness- and structure-wise, a continuous transition threshold function F(r,σ,k,h)=exp(−max(r 2 −k·σ 2 ,0)/h 2 ) is used, where k and h represent parameters. 
         [0086]    Vectorized version of combined filtering can be defined as follows: 
         [0000]    
       
         
           
             
               
                 
                   
                     
                       U 
                       ^ 
                     
                      
                     
                       ( 
                       
                         
                           p 
                           m 
                         
                         , 
                         f 
                       
                       ) 
                     
                   
                   = 
                   
                     
                       ( 
                       
                         
                           ∑ 
                           
                             q 
                             ∈ 
                             
                               S 
                                
                               
                                 ( 
                                 
                                   
                                     p 
                                     m 
                                   
                                   , 
                                   t 
                                 
                                 ) 
                               
                             
                           
                         
                          
                         
                           
                             W 
                              
                             
                               ( 
                               
                                 
                                   p 
                                   m 
                                 
                                 , 
                                 q 
                                 , 
                                 f 
                               
                               ) 
                             
                           
                           · 
                           
                             U 
                              
                             
                               ( 
                               
                                 q 
                                 , 
                                 f 
                               
                               ) 
                             
                           
                         
                       
                       ) 
                     
                     · 
                     
                       / 
                       
                         
                           ∑ 
                           
                             q 
                             ∈ 
                             
                               S 
                                
                               
                                 ( 
                                 
                                   
                                     p 
                                     m 
                                   
                                   , 
                                   t 
                                 
                                 ) 
                               
                             
                           
                         
                          
                         
                           
                             W 
                              
                             
                               ( 
                               
                                 
                                   p 
                                   m 
                                 
                                 , 
                                 q 
                                 , 
                                 f 
                               
                               ) 
                             
                           
                           . 
                         
                       
                     
                   
                 
               
               
                 
                   ( 
                   6 
                   ) 
                 
               
             
           
         
       
     
         [0087]    Here, W(p m ,q,f)=V(B(p m ,q,f),n(p m ,q)) represents block of values of weight function of combined filter; and B(p m ,q,f) represents block of values of weight function of the bilateral filter, B(p m ,q,f)={b(p m +s,q+s),sεS(0, f)}, and both such values are calculated as b(p m +s,q+s)=G(s,a)·F((ū(p m +s)−ū(q+s)),σ(ū(p m +s)),k b ,h b ), where G(s,a) represents a Gaussian function; n(p m ,q)=F(∥Ū(p m , f)−Ū(q, f)∥ 2 ,σ(ū(p m )),k n ,h n ) represents weight of non-local means filter, which is the same for all pixels of the comparison block; σ(u(p m )) and σ(Ū(p m , f)) represent the value and the RMS error block of estimation image pixel noise, dependent on signal intensity; • represents an operator for bit-by-bit multiplication/division of blocks; and Û(p m , f) represents the result estimation of intensities of pixels within the current partition block. 
         [0088]    In the combined filtering process, described here (5), (6), it is of importance to properly select weight function V that can be determined by various methods. In terms of execution speed-effective combined filtering algorithm, simple yet quite qualitatively effective version of combined filter can be created by framing weight function as a preset proportion of brightness and structural weights: 
         [0000]        V ( b,n )=(1−λ)· b ( p,q )+λ· n ( p,q );  (7)
 
         [0089]    and for vectorized version of combined filter: 
         [0000]        V ( B,n )=(1−λ)· B ( p   m   ,q,f )+λ· n ( p   m   ,q ).  (8)
 
         [0090]    The selected value of λ should be close to maximum possible non-local means filter weight value, e.g. λ≈0.8. Selection of (7) or (8) weight function allows separating bilateral filter and non-local filter. If filter separation is used, combined filtration comprises sequential bilateral and non-local filtering and further combining obtained results. Separation of filters allows selecting more speed-effective values of filter parameters (e.g. reducing bilateral filter averaging area), as well as eliminating unnecessary bilateral filtering of pixels, associated with regions where partition blocks overlap. By the example of vectorized combined filter (6), when using weight function (8), separable combined filtering can be written down as follows: 
         [0000]      {circumflex over ( U )}( p   m   ,f )==(λ· Û   n ( p   m   ,f )+(1−λ)· Û   b ( p   m   ,f ))/(λ· c ( p   m )+(1−λ)· C ( p   m   ,f ))  (9)
 
         [0091]    Here, 
         [0000]    
       
         
           
             
               
                 
                   U 
                   ^ 
                 
                 n 
               
                
               
                 ( 
                 
                   
                     p 
                     m 
                   
                   , 
                   f 
                 
                 ) 
               
             
             = 
             
               
                 ∑ 
                 
                   q 
                   ∈ 
                   
                     S 
                      
                     
                       ( 
                       
                         
                           p 
                           m 
                         
                         , 
                         t 
                       
                       ) 
                     
                   
                 
               
                
               
                 
                   w 
                    
                   
                     ( 
                     
                       
                         p 
                         m 
                       
                       , 
                       q 
                     
                     ) 
                   
                 
                  
                 
                   U 
                    
                   
                     ( 
                     
                       q 
                       , 
                       f 
                     
                     ) 
                   
                 
               
             
           
         
       
     
         [0000]    represents the result of filtering by vectorized non-local means filter (without normalizing with respect to total weight); Û b (p m , f) represents an image pixel block, processed by bilateral filter (without normalizing with respect to total weight), therefore 
         [0000]    
       
         
           
             
               
                 
                   
                     U 
                     ^ 
                   
                   b 
                 
                  
                 
                   ( 
                   
                     
                       p 
                       m 
                     
                     , 
                     f 
                   
                   ) 
                 
               
               = 
               
                 { 
                 
                   
                     
                       u 
                       ^ 
                     
                      
                     
                       ( 
                       s 
                       ) 
                     
                   
                   , 
                   
                     
                       
                         s 
                         ∈ 
                         
                           S 
                            
                           
                             ( 
                             
                               
                                 p 
                                 m 
                               
                               , 
                               f 
                             
                             ) 
                           
                         
                       
                       | 
                       
                         
                           u 
                           ^ 
                         
                          
                         
                           ( 
                           s 
                           ) 
                         
                       
                     
                     = 
                     
                       
                         ∑ 
                         
                           q 
                           ∈ 
                           
                             S 
                              
                             
                               ( 
                               
                                 s 
                                 , 
                                 t 
                               
                               ) 
                             
                           
                         
                       
                        
                       
                         
                           b 
                            
                           
                             ( 
                             
                               s 
                               , 
                               q 
                             
                             ) 
                           
                         
                          
                         
                           u 
                            
                           
                             ( 
                             q 
                             ) 
                           
                         
                       
                     
                   
                 
                 } 
               
             
             ; 
           
         
       
       
         
           
             
               c 
                
               
                 ( 
                 
                   p 
                   m 
                 
                 ) 
               
             
             , 
             
               C 
                
               
                 ( 
                 
                   
                     p 
                     m 
                   
                   , 
                   f 
                 
                 ) 
               
             
           
         
       
     
         [0000]    represents total weight and total weight block of non-local means filter and bilateral filter, respectively. 
         [0092]    When separable bilateral filter is used, separable form of combined filtering (9) also allows increasing filtering speed, making it considerably closer to execution speed of vectorized version of non-local means filter. 
         [0093]    Effectiveness of combined filtering method (6) with type (8) weight is demonstrated in  FIG. 5  (A, B, C, D, E), where:  5 A is a fragment of radiological test pattern;  5 B is an image that has been processed by vectorized version (4) of standard non-local means filter;  5 C is an image that has been processed by combined filter;  5 D is a difference between  5 A and  5 B; and  5 E is a difference between  5 A and  5 C. Artifacts, typical of standard non-local means filtering method, well visible along object borders (figure strips)), turn out to be considerably reduced as a result of application of combined filter. 
         [0094]    Stage 3. Low Frequency Noise Suppression. 
         [0095]    In smooth region of an image (i.e. areas of the same brightness with homoscedasticity manifesting itself as noise), bilateral filter, as well as non-local means filter, shows characteristics of an elementary linear neighborhood averaging filter (i.e. it operates as moving average). Due to using finite averaging region, part of spectral noise components, associated with low frequencies, will not be suppressed and will become well visible in smooth areas of an image. 
         [0096]    The present invention uses image pyramid with bilateral linear regression filtering for residual LF noise suppression [Burt et. al.,  The Laplacian pyramid as a compact image code, IEEE Trans. Commun ., vol. COM-31, pp. 532-540, April 1983; Tomasi et. al.,  Bilateral Filtering for Gray and Color Images, In Proc.  6 th Int. Conf Computer Vision , New Delhi, India, 1998, pp. 839-846]. Image pyramid is constructed based on the following formulae: 
         [0000]    
       
         
           
             
               
                 
                   { 
                   
                     
                       
                         
                           
                             
                               
                                 A 
                                 0 
                               
                                
                               
                                 ( 
                                 p 
                                 ) 
                               
                             
                             = 
                             
                               I 
                                
                               
                                 ( 
                                 p 
                                 ) 
                               
                             
                           
                           , 
                         
                       
                     
                     
                       
                         
                           
                             
                               
                                 A 
                                 k 
                               
                                
                               
                                 ( 
                                 p 
                                 ) 
                               
                             
                             = 
                             
                               Reduce 
                                
                               
                                 ( 
                                 
                                   
                                     A 
                                     
                                       k 
                                       - 
                                       1 
                                     
                                   
                                    
                                   
                                     ( 
                                     p 
                                     ) 
                                   
                                 
                                 ) 
                               
                             
                           
                           , 
                         
                       
                     
                     
                       
                         
                           
                             
                               
                                 H 
                                 
                                   k 
                                   - 
                                   1 
                                 
                               
                                
                               
                                 ( 
                                 p 
                                 ) 
                               
                             
                             = 
                             
                               
                                 
                                   A 
                                   
                                     k 
                                     - 
                                     1 
                                   
                                 
                                  
                                 
                                   ( 
                                   p 
                                   ) 
                                 
                               
                               - 
                               
                                 Expand 
                                  
                                 
                                   ( 
                                   
                                     
                                       A 
                                       k 
                                     
                                      
                                     
                                       ( 
                                       p 
                                       ) 
                                     
                                   
                                   ) 
                                 
                               
                             
                           
                           , 
                           
                             k 
                             = 
                             
                               
                                 
                                   1 
                                   , 
                                   n 
                                 
                                 _ 
                               
                               .. 
                             
                           
                         
                       
                     
                   
                 
               
               
                 
                   ( 
                   10 
                   ) 
                 
               
             
           
         
       
     
         [0097]    Here, I(p) represents an original image; A k (p),H k (p) represents low-frequency (approximating) and high-frequency (detailing) images of level k of the image pyramid; Reduce represents an operator for filtering and double reduction of image sampling rate (decimation), i.e. Reduce(A)=(A*L r )↓2, where L r  represents linear LF filter, and ↓2 represents decimation; Expand represents an operator for image downsampling and filtering, i.e. Expand(A)=(A↑2)*L e , where ↑2 represents downsampling, and L e  represents linear low-frequency filter. 
         [0098]    Noise variance distribution is estimated based on low-frequency components of the pyramid (10). To achieve that, white Gaussian noise (AWGN) with standard deviation σ is supplied to the combined filter input. Based on the above-mentioned equivalence of combined filter and neighborhood averaging filter in smooth image areas, dependence of LF images RMS error value on pyramid level number is calculated: 
         [0000]    
       
         
           
             
               
                 
                   { 
                   
                     
                       
                         
                           
                             
                               σ 
                               0 
                             
                             = 
                             
                               σ 
                               / 
                               
                                 ( 
                                 
                                   
                                     2 
                                     · 
                                     t 
                                   
                                   + 
                                   1 
                                 
                                 ) 
                               
                             
                           
                           , 
                         
                       
                     
                     
                       
                         
                           
                             
                               σ 
                               k 
                             
                             = 
                             
                               
                                 σ 
                                 0 
                               
                               · 
                               
                                 q 
                                  
                                 
                                   ( 
                                   
                                     k 
                                     , 
                                     
                                       L 
                                       r 
                                     
                                   
                                   ) 
                                 
                               
                             
                           
                           , 
                           
                             k 
                             = 
                             
                               
                                 1 
                                 , 
                                 n 
                               
                               _ 
                             
                           
                           , 
                         
                       
                     
                   
                 
               
               
                 
                   ( 
                   11 
                   ) 
                 
               
             
           
         
       
     
         [0099]    where σ 0  represents RMS error level of noise in an image, subjected to combined filtering; σ k  represents RMS error level of image noise A k (p) in pyramid (10); and q(k, L r ) represents reduction rate of RMS error of noise in images A k (p). 
         [0100]    To suppress residual LF noise, images A k (p) of pyramid (10) are processed by bilateral filtering. To reduce staircasing artifacts [Buades et. al.,  The staircasing effect in neighborhood filters and its solution, IEEE Transactions on Image Processing , Vol. 15 (6), pp. 1499-1505, 2006], bilateral linear regression filter is applied: 
         [0000]        â   k ( i,j )= a·i+b·j+c.   (12)
 
         [0101]    Parameters (a,b,c) are selected based on functional minimization requirement 
         [0000]    
       
         
           
             
               
                 min 
                 
                   a 
                   , 
                   b 
                   , 
                   c 
                 
               
                
               
                 ( 
                 
                   
                     ∑ 
                     
                       
                         ( 
                         
                           x 
                           , 
                           y 
                         
                         ) 
                       
                       ∈ 
                       
                         S 
                          
                         
                           ( 
                           
                             i 
                             , 
                             j 
                             , 
                             t 
                           
                           ) 
                         
                       
                     
                   
                    
                   
                     
                       w 
                        
                       
                         ( 
                         
                           x 
                           , 
                           y 
                         
                         ) 
                       
                     
                      
                     
                       
                         ( 
                         
                           
                             
                               a 
                               k 
                             
                              
                             
                               ( 
                               
                                 x 
                                 , 
                                 y 
                               
                               ) 
                             
                           
                           - 
                           
                             a 
                             · 
                             i 
                           
                           - 
                           
                             b 
                             · 
                             j 
                           
                           - 
                           c 
                         
                         ) 
                       
                       2 
                     
                   
                 
                 ) 
               
             
             ; 
           
         
       
     
         [0000]    S(i, j,t) represents an area of radius t, centered at the current pixel a k (i, j) of the image with number k of pyramid (10); w(x, y)=G(x−i, y−j,a)·F((a k (i, j)−a k (x, y)),σ k (i, j)) represents averaging weight, where spatial proximity of pixels is defined by Gaussian function G(x−i, y−j,a)=exp(−((x−i) 2 +(y−j) 2 )/a 2 ), and brightness proximity of pixels is estimated by function F((a k (i, j)−a k (x, y),σ k (i, j))=exp(−(a k (i, j)−a k (x, y)) 2 /σ k   2 (i, j)); σ k (i, j) represents signal-dependent RMS error of approximating pyramid image, defined by noise map and dependence of type (11); and â k (i, j) represents the result estimation of intensity of the current image pixel. To speed up filtering process, separable filter version (12) is used. 
         [0102]    Filtered approximating components of the pyramid are used to restore the image with reduced level of residual low-frequency noise 
         [0000]      {circumflex over ( A )} k-1 ( p )={circumflex over ( H )} k-1 ( p )+Expand({circumflex over ( A )} k ( p )), k =    n, 1 .  (13)
 
         [0103]    Here, symbols Â, Ĥ refer to the fact, that construction of pyramid (10) involved filtering of LF images by filter (12). 
         [0104]    Stage 4. Over-Smoothing Correction. 
         [0105]    In the present invention, over-smoothing of useful high-frequency information is corrected by blending filtered image and its original version together. Blending weights are calculated based on pixel and structural similarities of filtered image and original image. To determine structural similarity, Euclidean distance between pixel blocks of the images is calculated; to determine pixel similarity (intensity similarity), the distance between brightness levels of pixels of the relevant images is calculated. 
         [0106]    Let û(p) and u(p) represent intensity values of filtered and original image in pixel p, respectively. In this case, pixel similarity of frames is calculated based on the following threshold function: 
         [0000]        w   i ( p )=exp(−max(({circumflex over ( u )}( p )− u ( p )) 2   −k   i ·σ 2 ( u ( p )),0)/ h   i   2 ),  (14)
 
         [0107]    where k i , h i  represent parameters that define characteristic shape of distribution; and σ 2 (u(p)) represents noise variance value in pixel p. In turn, structural similarity of filtered image and original image is determined based on the following formula: 
         [0000]        w   s ( p )=exp(max(−(∥{circumflex over ( U )}( p )− U ( p )∥ 2   2   −k   s ·σ 2 ( u ( p )),0)/ h   s   2 ,  (15)
 
         [0108]    where k s , T s  represent parameters, defining shape of threshold function; and ∥Û(p)−U (p)∥ 2   2  represents squared Euclidian distance between pixel blocks of the respective images. 
         [0109]    Calculated pixel and structural similarities of the images are combined to determine the final measure of pixel-structural similarity in the following way: 
         [0000]        w ( p )= w   i ( p )· w   i ( p )+(1− w   i ( p ))· w   s ( p ).  (16)
 
         [0110]    In this case, over-smoothing correction involves combining filtered image and original noisy image based on the following rule: 
         [0000]        û   c ( p )= w ( p )·{circumflex over ( u )}( p )+(1− w ( p ))· u ( p ).  (17)
 
         [0111]    Setting similarity measure as (16), along with correction based on formula (17), allows filtering individual abrupt image changes (i.e. outliers). 
         [0112]    To estimate image noise variance, firstly, estimation image and noise image are obtained by applying to original image I(x, y) the following 3×3 low-frequency linear binomial filter H 1 =[1, 2, 1]/4, H=H 1   T ·H. It allows creating smoothed image I e (p)=I*H, where * represents convolution. After that, noise image N e (p)=I(p)−I e (p) is calculated. 
         [0113]    By removing borders from variance calculations based on noise image N e (p), images of positive and negative changes are formed: 
         [0000]    
       
         
           
             
               
                 N 
                 e 
                 - 
               
                
               
                 ( 
                 p 
                 ) 
               
             
             = 
             
               { 
               
                 
                   
                     
                       
                         
                           1 
                           , 
                         
                       
                       
                         
                           
                             
                               N 
                               e 
                             
                              
                             
                               ( 
                               p 
                               ) 
                             
                           
                           &lt; 
                           0 
                         
                       
                     
                     
                       
                         
                           0 
                           , 
                         
                       
                       
                         
                           
                             
                               
                                 N 
                                 e 
                               
                                
                               
                                 ( 
                                 p 
                                 ) 
                               
                             
                             ≥ 
                             0 
                           
                           , 
                         
                       
                     
                   
                    
                   
                     
 
                   
                    
                   
                     
                       N 
                       e 
                       + 
                     
                      
                     
                       ( 
                       p 
                       ) 
                     
                   
                 
                 = 
                 
                   { 
                   
                     
                       
                         
                           1 
                           , 
                         
                       
                       
                         
                           
                             
                               N 
                               e 
                             
                              
                             
                               ( 
                               p 
                               ) 
                             
                           
                           &gt; 
                           0 
                         
                       
                     
                     
                       
                         
                           0 
                           , 
                         
                       
                       
                         
                           
                             
                               N 
                               e 
                             
                              
                             
                               ( 
                               p 
                               ) 
                             
                           
                           ≤ 
                           0. 
                         
                       
                     
                   
                 
               
             
           
         
       
     
         [0114]    To separate large areas, associated with borders, these binary images are subjected to step-by-step morphological erosion and dilatation, using masks of 2×2 size: 
         [0000]        B   e   −1 ( p )=dilate [2×2] (erode [2×2 ]( N   e   − ( p ))), B   e   + ( p )=dilate [2×2 ](erode [2×2 ]( N   e   + ( p ))). 
         [0000]    Then these images are combined to create a map of the original image borders E(p)=B e   − (p)∩B e   + (p). 
         [0115]    To calculate interval estimations of noise variance, minimum I min  and maximum I max  of image intensity I e (p) are determined, increment h M  is selected, and intensity range is partitioned into intervals M i  in increments of h M  (e.g. h M =32). For each image pixel I e (p), an interval such pixel&#39;s value pertains to is identified, and the relevant pixel value of the image N e (p) is used to calculate noise variance estimation σ 2 (i) within that interval M i , which involves eliminating pixels with border map value E(p)=1. In calculating interval estimation of noise variance, unbiased variance estimator is used: 
         [0000]    
       
         
           
             
               
                 
                   σ 
                   2 
                 
                  
                 
                   ( 
                   i 
                   ) 
                 
               
               = 
               
                 
                   ∑ 
                   
                     j 
                     = 
                     1 
                   
                   
                     n 
                     i 
                   
                 
                  
                 
                   
                     
                       ( 
                       
                         
                           
                             N 
                             e 
                             i 
                           
                            
                           
                             ( 
                             j 
                             ) 
                           
                         
                         - 
                         
                           
                             N 
                             _ 
                           
                           e 
                           i 
                         
                       
                       ) 
                     
                     2 
                   
                   / 
                   
                     ( 
                     
                       
                         n 
                         i 
                       
                       - 
                       1 
                     
                     ) 
                   
                 
               
             
             , 
           
         
       
     
         [0116]    where N e   i (j) represents the value of pixel of noise image N e (p) from interval M i , n i  represents total amount of accumulated values of noise image pixels within interval M i , and  N   i   e  represents average value of noise pixels within interval M i . 
         [0117]    The obtained interval estimations of variance σ 2  (i) are refined in such a way as to iteratively eliminate in each interval M i  the values of noise pixels that exceed in absolute magnitude a threshold, equal to three standard noise deviations, which is followed by conversion of noise variance estimation in that interval: 
         [0000]        N   e   i   [k+ 1]={ N   e ( p )|( N   e ( p )ε N   e   i   [k])∩(   N   e ( p )≦3σ( i ))},
 
         [0118]    where N e   i [k+1] represents a set of values of noise pixels within interval M i  with iteration k. Further calculations use only intervals, which have accumulated sufficient number of values of noise pixels (e.g. at least 500). Besides, excluded from consideration are intervals with average value  N   i   e  significantly different from zero (deviating from zero by more than a half of interval grid increment h M ), as such intervals are most probably dominated by residual values of non-noise pixels. 
         [0119]    To estimate relationship of noise variance and intensity, interpolating tabulated function is created based on noise variance estimations obtained. Such tabulated function is formed based on robust local linear approximation of interval estimations of noise variance. For that, on interval grid M i , increment h I  and approximation radius r l  are selected (they can be made dependent on the number of points n i  within intervals M i ). Values of variance-intensity tabulated function, developed at the previous stage, are approximated based on the following rule: 
         [0000]      {circumflex over (σ)} 2 ( k·h   I )= a·m   k   +b,  
 
         [0120]    where k represents the number of approximation node; m k  represents the centre of interval M i (k·h I ); and parameters a,b are calculated based on sum of absolute deviations [Press et. al., Numerical Recipes in C: The Art of Scientific Computing, Second Edition. New York: Cambridge University Press, 1999] 
         [0000]    
       
         
           
             
               
                 ∑ 
                 
                   j 
                   = 
                   
                     
                       kh 
                       I 
                     
                     - 
                     
                       r 
                       I 
                     
                   
                 
                 
                   
                     kh 
                     I 
                   
                   + 
                   
                     r 
                     I 
                   
                 
               
                
               
                  
                 
                   
                     
                       
                         σ 
                         ^ 
                       
                       2 
                     
                      
                     
                       ( 
                       j 
                       ) 
                     
                   
                   - 
                   
                     a 
                     · 
                     
                       m 
                       k 
                     
                   
                   - 
                   b 
                 
                  
               
             
             -&gt; 
             
               
                 min 
                 
                   a 
                   , 
                   b 
                 
               
               . 
             
           
         
       
     
         [0121]    The resultant list of values {M i (k·h I ),{circumflex over (σ)} 2 (k·h I )} is interpolated over the entire range of intensities [I min ,I msc ], thus forming interpolating table of the desired dependence of noise variance on signal intensity. 
         [0122]    Further on, based on the estimation image and the interpolating table, a noise map is developed, i.e. an image where each pixel estimates noise variance in associated pixel of the original image. 
         [0123]    To ensure proper balance between filtering quality and filtering speed, combined bilateral non-local filtering is performed using formula (9), where λ=0.85, averaging region radius t=5, pixel block radius f=4, and blocks are assumed to overlap by 1 pixel. To determine similarity of pixels from the point of view of brightness and structure, a continuous transition threshold function F(r,σ,k,h)=exp(−max(r 2 −k·σ 2 ,0)/h 2 ) is used, where k and h represent parameters. Parameter value for brightness similarity is assumed as k b =9, and parameter value for structural similarity is assumed as k n =2.1. Values of parameter h b  are selected in such a way as to achieve the pre-set rate of decreasing of function F, which is achieved by setting function value F(r,σ,k,h), e.g. as F=0.1, at preset threshold value level kσ 2 , e.g. as kσ 2 =10·σ 2 , which results in parameter value h b =9. In calculation of weight function, geometric distance is not taken in consideration. In similar way, parameter h n  value for non-local means filter weight function is identified: at level kσ 2 =3·σ 2 , function value is assumed as F=0.5, which results in parameter value h n =1. Estimation image, used as a basis for calculation of averaging weights, is determined by filtering the original image with the following LF linear binomial filter of H 1 =[1, 2, 1]/4, H=H 1   T ˜H size. The value of parameter σ is taken from the noise map, developed at noise estimation stage. 
         [0124]    At residual LF noise filtering stage, image pyramid (10) with number of levels n=4 is used, wherein low frequency filter within the Reduce operator is assumed as L r =[1, 1]2, L r =L r   T ·L r . With this kind of filter, noise RMS error level with low-frequency components of the image pyramid decreases by law of q(k, L r )=2 −k , k=  1,n . In the course of construction of the pyramid HF components after downsampling, LF filter L e =[1,4,6,4,1]/16. L e =L e   T ·L e  is used in Expand operator. For LF filtering, bilateral filter is used, where geometric distance is not taken into account; filter radius is assumed as 5. Weight function parameters coincide with those of weight function of bilateral filter in the combined filter. In calculation of weight function, geometric distance is not taken in consideration. 
         [0125]    In over-smoothing correction, to calculate structural similarity, blocks of radius  5  are used. Parameters, defining shape of weight functions, coincide with weight function parameters of the combined filter. In calculation of weight function of intensity-based pixel similarity, geometric distance is not taken in consideration. 
       REFERENCES 
       [0000]    
       
         1. Bilcu R. C., Vehvilainen M. Fast non-local means for image de-noising. In Digital Photography III Proc. of SPIE IS&amp;T, Russel A. Martin, Jeffrey M. DiCarlo, and Nitin Sampat, Eds., February 2007, vol. 6502. 
         2. Bosco A., Battiato S., Bruna A. R., Rizzo R. Noise Reduction for Cfa Image Sensors Exploiting Hvs Behavior. Sensors 9(3), 1692-1713, 2009. 
         3. Buades A., Coll B., Morel J.-M. A review of image denoising algorithms, with a new one. Multiscale Modeling and Simulation (SIAM interdisciplinary journal), Vol. 4 (2), pp. 490-530, 2005. 
         4. Buades A., Coll B., Morel J.-M. Image Denoising by Non-Local Averaging. Proc. of IEEE Int. Conf. on Acoustics, Speech, and Signal Processing, vol. 2, pp. 25-28, March 2005. 
         5. Buades A., Coll B., Morel J.-M. On image denoising methods. SIAM Multiscale Modeling and Simulation, 4(2):490-530, 2005. 
         6. Buades A., Coll B., Morel J.-M. The staircasing effect in neighborhood filters and its solution. IEEE Transactions on Image Processing, Vol. 15 (6), pp. 1499-1505, 2006. 
         7. Burt P., Adelson E. The Laplacian pyramid as a compact image code. IEEE Trans. Commun., vol. COM-31, pp. 532-540, April 1983. 
         8. Candes E. J., Donoho D. Curvelets, multiresolution representation, and scaling laws. In SPIE conference on Signal and Image Processing: Wavelet Applications in Signal and Image Processing VIII, San Diego, 2000. 
         9. Dauwe A., Goossens B., Luong H., Philips W. A Fast Non-Local Image Denoising Algorithm. In Proc. SPIE Electronic Imaging, vol. 6812, San Jose, USA, January 2008. 
         10. Foi A., Trimeche M., Katkovnik V., Egiazarian K. Practical poissonian-gaussian noise modeling and fitting for single-image raw-data. Image Processing, IEEE Transactions on 17, 1737-1754 (October 2008). 
         11. Foi A. Practical Denoising of Clipped or Overexposed Noisy Images. Proc. 16th European Signal Process. Conf., EUSIPCO 2008, Lausanne, Switzerland, August 2008. 
         12. Forstner W. Image Preprocessing for Feature Extraction in Digital Intensity, Color and Range images. In Springer Lecture Notes on Earth Sciences, 1998. 
         13. Gino M. Noise, Noise, Noise. http://www.astrophys-assist.com/educate/noise/noise.htm 
         14. Hensel M., Lundt B., Pralow T., Grigat R.-R. Robust and Fast Estimation of Signal-Dependent Noise in Medical X-Ray Image Sequences. In Handels, H., et al., ed.: Bildverarbeitung fur die Medizin 2006: Algorithmen, Systeme, Anwendun-gen. Springer (2006) 46-50. 
         15. Kervrann C., Bourlanger J. Optimal spatial adaptation for patch-based image denoising. IEEE Transactions on Image Processing, vol. 15, no. 10, October 2006. 
         16. Kingsbury N. G. The dual-tree complex wavelet transform: A new technique for shift invariance and directional filters. In Proc. 8th IEEE DSP Workshop, Utah, Aug. 9-12, 1998, paper no. 86. 
         17. Liu C., Szeliski R., Kang S. B., Zitnick C. L., Freeman W. T. Automatic estimation and removal of noise from a single image. IEEE Transactions on Pattern Analysis and Machine Intelligence 30(2), 299-314, 2008. 
         18. Mahmoudi M., Sapiro G. Fast image and video denoising via nonlocal means of similar neighborhoods. IEEE Signal Processing Letters, 12(12):839-842, 2005. 
         19. Peter J., Govindan V. K., Mathew A. T. Robust estimation approach for nonlocal-means denoising based on structurally similar patches. Int. J. Open Problems Compt. Math., vol. 2, no. 2, June 2009. 
         20. Press W., Teukolsky S., Vetterling W., Flannery B. Numerical Recipes in C: The Art of Scientific Computing, Second Edition. New York: Cambridge University Press, 1992. 
         21. Rudin L., Osher S., Fatemi E. Nonlinear total variation based noise removal algorithms. Physica D, 60:259-268, 1992. 
         22. Salmeri M., Mencattini A., Rabottino G., Lojacono R. Signal-dependent Noise Characterization for Mammographic Images Denoising. IMEKO TC4 Symposium (IMEKOTC4 &#39;08), Firenze, Italy, September 2008. 
         23. Selesnick I. W. The double-density dual-tree discrete wavelet transform. IEEE Trans. Signal Processing, vol. 52, no. 5, pp. 1304-1314, May 2004. 
         24. Starck J., Murtagh F., Bijaoui A. Image Processing and Data Analysis: The Multiscale Approach. Cambridge University Press, 1998. 
         25. Strang G., Nguyen T. Wavelets and Filter Banks. Wellesley-Cambridge Press. 
         26. Takeda H., Farsiu S., Milanfar P. Kernel regression for image processing and reconstruction. IEEE Transactions on Image Processing, vol. 16, no. 2, pp. 349-366, February 2007. 
         27. Tomasi C., Manduchi R. Bilateral Filtering for Gray and Color Images. In Proc. 6th Int. Conf. Computer Vision, New Delhi, India, 1998, pp. 839-846. 
         28. Vetterli M., Kovacevic J. Wavelets and Subband Coding. Prentice Hall, Englewood Cliffs, N.J., 1995. 
         29. Waegli B. Investigations into the Noise Characteristics of Digitized Aerial Images. In: Int. Arch. For Photogr. And Remote Sensing, Vol. 32-2, pages 341-348, 1998. 
         30. Wang J., Guo Y., Ying Y., Liu Y., Peng Q. Fast non-local algorithm for image denoising. In Proc. Of IEEE International Conference on Image Processing (ICIP), 2006, pp. 1429-1432. 
         31. Weickert J. Anisotropic Diffusion in Image Processing. Tuebner Stuttgart, 1998. 
         32. Yane B., Digital image processing. Moscow, Technosphera, 2007, p. 583