Abstract:
The present invention uses five regions of interest, 4 corners and 1 center to do Fourier Transform analysis to mark rough location of the streaks, if there is any. It sets the small window around the center in each Fourier Transform ROI to zero to mask the cluster of bright points caused by the lens rolloff or other noises. It use iterative linear regression to remove the random outliers and to search the best possible linear pattern. If it finds any streaks, it calculates the angle and converts it to the streak angle in spatial domain. Once a streak is detected, it can be removed by applying an inverse Fourier Transform on the processed magnitude and original phase Fourier Transform images.

Description:
FIELD OF THE INVENTION 
     The invention relates generally to the field of detecting and removing noise, such as streaks, from images and the like. More specifically, the invention relates to applying a frequency transform, such as a Fourier transform, to an image for detecting streaks in the image which may be quantified and removed. The present invention also includes an improved method of detecting streaks in images, especially flat-field images or images having known transform characteristics. 
     BACKGROUND OF THE INVENTION 
     Prior art algorithms used an edge detection algorithm in the spatial domain to identify discontinuities in the image data and measure the change in intensity taking place over a number of pixels. It is difficult to set the threshold to detect bright pixels which may contribute to an edge due to the noise in the image. It is also difficult to set a region of interest and the criteria to calculate the defect density; susceptible to false detection due to lens roll off or edge effect. 
     U.S. patent application Ser. No. 11/180,816, filed Jul. 13, 2005, entitled “The Use Of Frequency Transforms In The Analysis Of Image Sensors”, by Greg L. Archer, et al. uses a single linear regression algorithm in the frequency domain to identify a linear pattern in the center of the image. 
     Although the currently known and used methods for detecting and removing streaks is satisfactory, improvements are always desirable. The present invention provides such an improvement. 
     SUMMARY OF THE INVENTION 
     The present invention is directed to overcoming one or more of the problems set forth above. Briefly summarized, according to one aspect of the present invention, the invention uses five regions of interest, 4 corners and 1 center to do Fourier Transform analysis to mark rough location of the streaks, if there is any. It sets the small window around the center in each Fourier Transform ROI to zero to mask the cluster of bright points caused by the lens rolloff or other noises. It use iterative linear regression to remove the random outliers and to search the best possible linear pattern. If it finds any streaks, it calculates the angle and converts it to the streak angle in spatial domain. Once a streak is detected, it can be removed by applying an inverse Fourier Transform on the processed magnitude and original phase Fourier Transform images. 
     Advantageous Effect of the Invention 
     The present invention has the following advantage of applying a frequency transform to an image and then applying an iterative regression technique for improved streak noise detection and removal. 
     These and other aspects, objects, features and advantages of the present invention will be more clearly understood and appreciated from a review of the following detailed description of the preferred embodiments and appended claims, and by reference to the accompanying drawings. 
    
    
     
       BRIEF DESCRIPTION OF THE DRAWINGS 
         FIG. 1  is a typical image containing noise; 
         FIG. 2  is a typical image with a streak type noise; 
         FIG. 3  is the resulting magnitude image of  FIG. 2  after applying the discrete Fourier transform; 
         FIGS. 4A-D  are flow charts of the software program of the present invention; 
         FIGS. 5A-C  are illustrations of the regression technique of the present invention; 
         FIG. 6  is a digital camera for illustrating a typical commercial embodiment for which the present invention is used. 
     
    
    
     DETAILED DESCRIPTION OF THE INVENTION 
     In the following description, the present invention will be described in the preferred embodiment as a software program. Those skilled in the art will readily recognize that the equivalent of such software may also be constructed in hardware. 
     It is also noted that the present invention may be used for image evaluation, image sensor testing, and image manipulation in addition to detecting and removing streaks. 
     Referring to  FIG. 1 , there is shown a typical image  10  having noise  20  therein. The noise  20  may be perceptible to the human eye or substantially non-perceptible to the human eye in which case of non-perceptibility the image will appear grainy or the like. The present invention detects and removes streaks from the image as described herein below. 
     An example of an image  30  with a streak type noise  40  is shown in  FIG. 2 . The image  30  illustrates a streaking noise  40 . Applying a discrete Fourier transform to the image of  FIG. 2  results in the magnitude image  50  shown in  FIG. 3  in the frequency domain. In this magnitude image, a line pattern  60  is shown positioned across the center which corresponds to the noise  40  in the spatial domain. The image is separated into four quadrants for facilitating the operation of the software program of the present invention. 
       FIGS. 4A-4D  are the flowcharts of the software program of the present invention. Referring to  FIG. 4A , the program is started S 2  and the resolution of the digital image  30  is obtained for determining whether the image size or resolution is a power of two S 4 . If it is greater than 1024×1024, the whole image is separated into five windows S 6  and each window size is 1024×1024. The five windows are the top left, bottom left, top right, bottom right and the center. Otherwise, it uses only one center window S 8  and its window size is chosen to the nearest size of a power of two. If original size is not of a power of 2, zero padding is required. All the windows need to be processed S 10  as described in the following paragraph. After one window is processed, check if all of them are done S 12 . If not, select next window S 14  and process again S 10  until all windows are processed. Then it outputs results S 16  and ends the program S 18 . 
     Referring to  FIG. 4B , the process (referred to as S 10 ) block starts S 20  with applying a Fourier transform on the image in selected window S 22 . Once the magnitude image  50  of the transform is obtained, the center region in the magnitude image  50  is set to a predefined value, preferably zero, S 24 . The purpose of this operation is to remove the noise points around the center. The size of the region is determined by the noise level, preferably 1/16 th  of the size of 4 regions which will be discussed in the next paragraph. The region can be a circle, a square or an eclipse or a rectangle or any other shapes. A logrithm operation is then applied S 26  to the whole magnitude image  50  to increase the contrast of the image. A detailed description on the logarithm operation will be described in detail herein below. 
     There are four regions (as referred to in  FIG. 3 ) around the center point, top right, top left, bottom right and bottom left. Since the top right and bottom left are mirror images, the same as top left and bottom right, preferably only two regions need to be processed, for example the top left and top right regions. First the top left region is selected S 28 . Then all the pixel values are compared against a predefined threshold value Td 1 , preferably the mean plus 3 times of standard deviation of the all pixel values inside the region, S 30 . If any pixel value is above this threshold, it is marked as a bright point. A linear regression process is applied on all marked bright points S 32  (as will be discussed in detail in  FIG. 4C ). After that, check if all regions have been selected S 34 . If no, select the top right region S 36  and repeat the procedures of S 30  and S 32 . After all regions have been processed, it reaches the end of process S 38 . 
     Referring to  FIG. 4C , the regression block starts S 40  and its regression counter resets to 0. The program first checks if the total marked points are above a predefined threshold Td 2 , preferably 5 and above, S 42 . If the total bright points are fewer than 5, no streaks are found S 44  and it reaches the end of regression S 46 . If there are 5 or more bright points, a first linear regression is applied to the marked points S 48  and regression count increases by 1. Then the R-Square value of the best fit line from the regression is calculated. If it is above another predefined threshold Td 3 , preferably 0.8, S 50 , a streak is found and its angle is calculated S 52  and it reaches the end of regression S 54 . A detailed description on how to calculate the streak angle is given later ( FIGS. 5A-C ). 
     Referring to  FIGS. 4C and 5A , if R-Squared of the fit line  70  is less than the Td 3 , all the distances between the marked points ( 1 - 11  in  FIG. 5A ) to the fit line  70  are calculated S 56 . The first point is selected S 58 . If its distance to the fit line is greater than a threshold value Td 4 , preferably the mean of the all distances, S 60 , the point is considered an outlier and is removed from the original group of all marked points S 62 . On the other hand, if the distance is within the threshold Td 4 , the point is kept in the group S 64 . After that, the program checks if all points have been processed S 66 . If not, the next point is selected and the whole comparison process is repeated S 68  until all the points are processed and the new subset of points are determined S 70 . This concludes one cycle of regression. If the number of regression is greater than the MaxIteration, preferably 5 iterations, S 72 , no streak is found S 74  and it ends the regression S 76 . If the number of regression is smaller than the MaxIteration, the program will use the new subset of points as the marked points (illustrated as points  1 ,  2 ,  3 ,  5 ,  6 ,  7  and  10  in  FIG. 5B ) and repeat the process from S 42 . The other points ( 4 ,  8 ,  9  and  11 ) are removed because their distance is greater than the threshold of Td 4 . 
     Referring to  FIG. 5B , another fit line  74  is calculated from the new subset of points. Assuming the R-Squared of the fit line  74  is less than Td 3 , then the distances from the remaining points to the fit line  74  are calculated. This distance is compared to the threshold Td 4 ′ (the mean of all the new distances, theoretically smaller than Td 4 ). For all distances smaller than Td 4 ′, these points are kept for the new subset of points (points  1 ,  2 ,  5 ,  6 , and  7 ) for further processing. If the R-Squared of the fit line  76  from the latest subset points ( 1 ,  2 ,  5 ,  6 , and  7 ) is greater than Td 3 , then the regression is finished and a streak is found. Referring to  FIG. 5C , the slope of the fit line  76  from the remaining points is calculated. The angle of the streak is calculated based on this slope  76 , i.e., perpendicular to the angle of the fit line  76 . 
     Referring to  FIG. 4D , another embodiment of this invention is to apply a fixed and predefined number of regression cycles, preferably 5, and then to evaluate the streak angle. It is noted for clarity of understanding that  FIG. 4D  is an alternative embodiment of  FIG. 4C . It starts S 80  and resets the iteration counter to 0. It uses the first set of marked points to compare a predefined threshold Td 2  S 82 . If the marked points are fewer than the Td 2 , no streaks are found S 84  and the regression session will end S 86 . If the points are more than the Td 2 , the first linear regression is applied S 88  and the iteration counter advances by 1. Then the distances from all the points to the fit line are calculated S 90 . After that, the first point is selected S 92  and its distance to the fit line is compared to another threshold Td 4  S 94 . If the distance is greater than Td 4 , the point is removed from the original group S 98 . If it is not, the point is kept in the ground S 96 . After that, the program checks if all points have been processed S 100 . If not, next point is selected S 102  and the whole comparison process is repeated S 94  until all the points are processed and the new subset of points are determined S 104 . This concludes one cycle of regression. If the number of regressions is smaller than the IterationNum (5) S 106 , the program will use the new subset of points as the marked points and repeat the process from S 82 . If, on the other hand, the number of regression equals IterationNum, the final R-Squared of the fit line is calculated S 108  and is compared to a predefined value Td 3 . If it is less than Td 3 , no streak is found S 110  and the regression ends S 114 . However, the R-Squared is greater than Td 3 , the streak angle is calculated S 112  followed by the end of the regression S 114 . 
     The Fourier transform used in the present invention S 22  is preferably applied to each row and then applied again vertically to each column. It is noted that Fourier transform is the preferred embodiment, but other transform methods may also be used. The two transforms together result in a magnitude image or plot. Typically, Fourier transforms that may be used are, but not limited to, fast Fourier transform (Fourier Transform) by Cooley and Tukey and discrete fast Fourier transforms by Danielson-Lanczos. The use of the Fourier transform provides the separation of the frequency content in the original image. Low frequency values represent little or no change in the image (overall shape), while high frequency values indicate rapid changes in the image over a short distance (details). Because there are discrete, equally spaced pixels in a digital image, it becomes efficient to use a Discrete Fourier Transform (DFT) version of the general Fourier transform. 
     The result of the DFT will be a magnitude image and a phase image. The magnitude image will provide equally spaced data representing the frequency domain. Different frequencies are represented at different distances from the origin. The value at the origin represents the DC component or average value of the original image data while values off of the origin represent different orientations in the original image. The pixel value or energy in the frequency domain indicates how much of that frequency and orientation is present in the original image. 
     The DFT calculation is shown as: 
     
       
         
           
             
               
                 
                   
                     
                       X 
                       ⁡ 
                       
                         ( 
                         k 
                         ) 
                       
                     
                     = 
                     
                       
                         1 
                         N 
                       
                       ⁢ 
                       
                         
                           ∑ 
                           
                             n 
                             = 
                             0 
                           
                           
                             N 
                             - 
                             1 
                           
                         
                         ⁢ 
                         
                           
                             x 
                             ⁡ 
                             
                               ( 
                               n 
                               ) 
                             
                           
                           ⁢ 
                           
                             ⅇ 
                             
                               
                                 - 
                                 j 
                               
                               ⁢ 
                               
                                 
                                   2 
                                   ⁢ 
                                   π 
                                   ⁢ 
                                   
                                       
                                   
                                   ⁢ 
                                   nk 
                                 
                                 N 
                               
                             
                           
                         
                       
                     
                   
                 
                 
                   
                     
                       
                         for 
                         ⁢ 
                         
                             
                         
                         ⁢ 
                         k 
                       
                       = 
                       0 
                     
                     , 
                     1 
                     , 
                     2 
                     , 
                     … 
                   
                 
               
             
             ⁢ 
             
                 
             
             , 
             
               N 
               - 
               1 
             
           
         
       
       
         
           
             N 
             = 
             
               # 
               ⁢ 
               
                   
               
               ⁢ 
               data 
               ⁢ 
               
                   
               
               ⁢ 
               samples 
             
           
         
       
       
         
           
             
               
                 
                   
                     x 
                     ⁡ 
                     
                       ( 
                       n 
                       ) 
                     
                   
                   = 
                   data 
                 
               
               
                 
                   
                     
                       for 
                       ⁢ 
                       
                           
                       
                       ⁢ 
                       n 
                     
                     = 
                     0 
                   
                   , 
                   1 
                   , 
                   2 
                   , 
                   … 
                   ⁢ 
                   
                       
                   
                   , 
                   
                     N 
                     - 
                     1 
                   
                 
               
             
           
         
       
       
         
           
             
               ⅇ 
               
                 
                   - 
                   j 
                 
                 ⁢ 
                 
                   
                     2 
                     ⁢ 
                     π 
                     ⁢ 
                     
                         
                     
                     ⁢ 
                     nk 
                   
                   N 
                 
               
             
             = 
             
               
                 cos 
                 ⁢ 
                 
                   
                     2 
                     ⁢ 
                     π 
                     ⁢ 
                     
                         
                     
                     ⁢ 
                     nk 
                   
                   N 
                 
               
               - 
               
                 j 
                 ⁢ 
                 
                     
                 
                 ⁢ 
                 sin 
                 ⁢ 
                 
                   
                     2 
                     ⁢ 
                     π 
                     ⁢ 
                     
                         
                     
                     ⁢ 
                     nk 
                   
                   N 
                 
               
             
           
         
       
     
     The exponent has an imaginary (j) term, making the result complex (real and imaginary values), even though the original image values were real. The magnitude and phase are calculated as follows: 
     
       
         
           
             magnitude 
             = 
             
               
                  
                 
                   X 
                   ⁡ 
                   
                     ( 
                     k 
                     ) 
                   
                 
                  
               
               = 
               
                 
                   
                     
                       
                         x 
                         r 
                         2 
                       
                       + 
                       
                         x 
                         i 
                         2 
                       
                     
                   
                   ⁢ 
                   
                       
                   
                   ⁢ 
                   phase 
                 
                 = 
                 
                   
                     tan 
                     
                       - 
                       1 
                     
                   
                   ⁡ 
                   
                     ( 
                     
                       
                         x 
                         i 
                       
                       
                         x 
                         r 
                       
                     
                     ) 
                   
                 
               
             
           
         
       
     
     An area image sensor with dimension M×N is illuminated with uniform light. The image data f(x,y) from the image sensor is an M×N array of voltages, one for each pixel from the image sensor. 
     The two-dimensional DFT of an image f(x,y) of size M×N is given by the equation 
               F   ⁡     (     u   ,   v     )       =       1   MN     ⁢       ∑     x   =   0       M   -   1       ⁢       ∑     y   =   0       N   -   1       ⁢       f   ⁡     (     x   ,   y     )       ⁢     ⅇ       -   j     ⁢           ⁢   2   ⁢           ⁢     π   ⁡     (       ux   /   M     +     vy   /   N       )                         
and the magnitude image equals the absolute value of F(u,v).
 
     In order to speed up the DFT calculation, a variant of DFT called fast Fourier transform (Fourier Transform) was developed in 1965. The requirement for utilizing the Fourier Transform algorithm is the x and y dimensions must each be a power of 2, such as 256×256, 512×512, 256×512, and 1024×1024. If the image size is not perfect power of 2, zero or average padding is needed. In this example, a 1024×1024 window region of interest is selected. 
     After the Fourier Transform of the image is obtained, the pixel in the frequency domain that represents the average value in the original image is set to zero. Then a logarithm transform (S 26 ) is applied to the magnitude image to boost the low magnitudes of peaks related to various frequencies. The transform is given by
 
 S ( u,v )=1 n (1 +|F ( u,v )|)
 
     In this example, the threshold value is the median value of S(u,v) in the quadrant. Then a linear regression of these thresholded values is performed based on the theory that a line should be present in the S(u, v) data which is perpendicular to the streak in the original captured image. The linear regression calculates a coefficient of correlation, R, which is compared with a pre-defined threshold value Td 3 . If R from the linear regression is greater than the predefined threshold value Td 3 , then the streak is detected. The fit line angle is calculated based on the equation angle=atan(slope) where the slope is calculated from the linear regression. 
     The streak angle is calculated (S 52  and S 112 ) due to the fact the streak angle is always perpendicular to the fit line angle. 
     Once a streak is detected based on the marked points at the end of regression, the points are preferably set to zero or alternatively substantially zero. Then the modified magnitude image in frequency domain along with the original phase image can be reversely transformed back to an original spatial image with the detected streak removed. 
     The digital image can be obtained by various methods. It can be captured by any imaging devices like a camera with an image sensor inside, an image scanner and etc. It can also come from a traditional film camera or film X-ray machine. After the image is formed on these traditional film medium, it can be digitized to a digital image. 
     Referring to  FIG. 6 , there is shown a digital camera  80  having the software program of the present invention installed in memory  90  and processed by a digital signal processor  100  for detecting and removing noise after image capture. This illustrates one of the above-described embodiments. 
     The invention has been described with reference to two preferred embodiments. However, it will be appreciated that variations and modifications can be effected by a person of ordinary skill in the art without departing from the scope of the invention. 
     PARTS LIST 
     
         
           10  Normal flat field image with random noise 
           20  Random noise 
           30  Flat field image with streak noise 
           40  Streak Noise 
           50  FFT Magnitude Image in Frequency Domain 
           60  A line pattern corresponding to the streak noise 
           70  The fit line of 1 st  regression 
           74  The fit line of 2 nd  regression 
           76  The final fit line 
           80  Digital camera 
           90  Memory 
           100  Digital signal processor 
         S 2 -S 114  Flowchart steps