Shot noise adaptive bilateral filter

Shot Noise Adaptive Bilateral Filter Harold Phelippeau(1),(2) , Hugues Talbot(1) , Mohamed Akil(1) ,Stefan Bara(2) (1)Université Paris-Est, Labinfo,ESIEE 93162 Noisy-le-Grand Cedex France {phelipph, talboth, akilm}@esiee.fr (2)NXP Semiconductors,2 Esplanade Anton Philips Campus Effiscience, Colombelles BP 2000 14906, Caen Cedex 9, France {harold.phelippeau, stefan.bara}@nxp.com Abstract 2. Bilateral filter characteristics The bilateral filter (BF) [1, 2, 3] is an important image denoising local filter. It reduces noise in images while preserving edges by means of nonlinear combination of local pixel values. Its formulation and implementation are both simple. However, the BF is not parameter-free. The set of the bilateral filter parameters has an important influence on its behavior and performance. They have to be chosen considering the end application. In the case of noise removal, the parameters have to be adapted to the noise level, while the bilateral filter adapts itself to the image details content. In this paper we propose a method to estimate the best bilateral filter intensity parameter set in the case of shot noise removal, the dominant form of natural noise in digital imaging. Index Terms— image sensors; Poisson noise model; noise filtering; noise reduction; noise estimation; photons flow; gaussian weight. 1. Introduction The bilateral filter is a non-linear filter well suited for denoising applications. It exhibits demonstrated effectiveness properties and its formulation simplicity contributes to its popularity. The Gaussian bilateral filter version, about which we focus this paper, has a set of parameters that have an important impact on filtering behavior and performance. We first recall the bilateral filter characteristics and the importance of the setting parameters. We then propose a method to automatically estimate the best intensity parameter set for shot noise removal, in the case of digital photography applications, by considering automatic gain control adjustment and the intensity mean and variance shot noise properties. We validate our assumptions with experiments on simulated sensor images. Finally we demonstrate the noise filtering adaptive advantage of the proposed method. (a) h=0 (b) h=10 (c) h=20 (d) h=40 Fig. 1. Bilateral filtering using a 5x5 square window as β and a variable h The bilateral filter replaces a pixel value in an image by a weighted mean of its neighbors considering both their geometric closeness and photometric similarities [1, 2, 3]. In this paper, we focus on the Gaussian bilateral filter because all the practical applications use this version [4]. The Gaussian bilateral filter is defined as follows: |u(y) − u(x)| 1 X − |x − y| exp − u(y) exp C(x) ρ2 h2 2 v(x) = 2 β (1) Where β represents the sliding window, y is a set of 2-D pixel positions in the sliding window, and x is the 2-D position of the centered pixel in the sliding window. u(x) is the intensity of the pixel at the x position in the original image, v(x) is the estimated pixel at the x position, ρ and h are respectively the standard deviation of the Gaussian distribution of the geometrical and the intensity weight. The filter behavior depends highly on these parameters setting. Parameter ρ can be chosen considering the size of the convolution kernel. Parameter h has to be chosen considering the level of filtering needed for the application. Indeed, the higher the standard deviation, the more the filter behaves like a low-pass Gaussian convolution. An illustration of the h influence is shown on Fig.1. For noise removal applications, the parameter h has to be chosen considering noise levels and statistics. Ideally, the best h should be the one that yields the lowest mean square error between the denoised image and the original version, however in real applications the original, noise-free image is never available. A possible solution consists of calibrating the best h for each noise level. Using this, finding the best h amounts to evaluating the noise level in the image. In the following section we present a novel method for shot noise level estimation. 3.3. Photons density from image statistics In the presence of gain manipulation only, we expect the noise variance to be proportional to average pixels values in nearconstant regions, due to its Poisson characteristics. Let consider I the image before AGC, and IG the adjusted image after AGC. Assuming G to be the gain factor, we can write: IG = G × I Considering a Poisson distribution process, we can assume the following relation where T − (I) represents regions of near constant values in the image I. E(T − (I)) = V ar(T − (I)) Best intensity parameter estimation based on image content 3. (2) (3) From image IG we can write: V ar(T − (IG )) = E(T − (IG )2 ) − (E(T − (IG )))2 (4) 3.1. Sensor noise assumption Our application is denoising in the context of digital photography. In this context, many sources are cause of noise generation in CCD and CMOS sensors. These can be categorized in four main types [5]: (1) The photon shot noise – associated with a random Poisson process governing the number of incident photons reaching a photosite; (2) the Photon Response Non Uniformity – caused by small sensitivity differences between photosites; (3) the dark current noise – produced by minority carriers thermally generated in the sensor wells; and (4) the read-out noise – resulting from thermal noise caused by MOSFET amplifiers. A sensor noise model complete in respect to these factors is proposed in [5]. According to this model, photon shot noise has, in standard illumination conditions, the most important influence on the output image. 3.2. Automatic gain control In a digital imaging system, an automatic gain control (AGC) is necessary to capture images in varying light conditions. It works by adjusting the average intensity of the output signal (see Fig.2). Assuming access to the gain factor and knowing the sensor resolution, it is possible to deduce the number of photons that have reached the sensor and then estimate the noise level. However, this information is not usually available to photographers, even using a RAW format. In the following sections we propose an estimation method of the AGC factor by using image statistics and the properties of a Poisson distribution process. Fig. 2. Digital conversion of a photon flow by a camera sensor, the AGC adjusts the average image output intensity Using (2) we can then derive: T − (IG ) = G × T − (I) (5) E(T − (IG )) = G × E(T − (I)) (6) and Replacing the expression (6) in (3) and (4), the gain is expressed: V ar(T − (IG )) (7) G= E(T − (IG )) The image mean E(IG ) and variance V ar(IG ) can be simultaneously calculated at each point in IG by using an efficient sliding-window algorithm [6]. We seek to establish a correlation between all E(IG ) and V ar(IG ) at each point where T − is non-zero. Since mean and variance are expected to be proportional, a robust linear, least-square correlation can be performed, the slope coefficient can then be interpreted as the gain factor transforming photon numbers to recorded pixel values. 3.4. Best h from photon density In order to estimate the gain, we need to select areas T − (IG ) in the noisy image IG that correspond to approximately constant values. These can be readily estimated from the image gradient T − (IG ) = εw (∇(Gσ ⋆ IG ) > ϑ), (8) where IG is the image, Gσ is a Gaussian convolution of variance σ, εw the morphological erosion on a window of size w × w ; ϑ is an intensity parameter. In our experiments, the gain could be accurately estimated by sampling V ar(IG ) and E(IG ) at 104 pixel locations. We chose σ = 2, ϑ = 2 and w = 11. As stated earlier, the bilateral filter smoothing properties vary with h. Intuitively we can assume that a low h should be preferable in low noise level, whereas a high h may be necessary in high-noise conditions. A calibration of the best h for each noise level is necessary to establish a correspondence between the gain factor estimation and the adapted h. This can be done using a digital imaging model or experimentally by calibrating sensors under different light illumination conditions. 4. Experiments and results In this section we seek to validate the assumptions of the previous section. 4.2. Photon number estimation We ran trials with a photon density per pixel ranging from 20 to 105 , i.e. from an extremely noisy to a very clean image, corresponding to over 10 f-stop, or alternatively from sensor sensitivity ranging from 6400 to 25 ASA. We estimated the number of photons Φ̄ with the following formula: P nbpix(I)E(I) (x,y)∈I I(x, y) Φ̄(I) = = , (9) G G where nbpix(I) is the total number of pixels in I. G is the estimated gain control, calculated with the method proposed in section 3.3. The correlation between simulated and estimated photon density is excellent, see Fig. 4 (R2 = 0.996). 8.5 ● 8.0 ● 7.5 ● 7.0 ● ● 6.5 ● ● 6.0 To do so we simulate image sensor acquisition via a regionalized cumulative spatial Poisson point process [7] using scene images as probability distribution functions. This simulates individual photons being recorded at each pixel location, as in a Monte-Carlo process. As scene images we used the Kodak PhotoCD database, extracted at the 512 × 768 resolution. These virtually noise-free images were scanned at the 3 full color samples per pixel from film original following the Kodak professional PhotoCD procedure. Figure 3 illustrates this simulation. estimated number of photons (log10 scale) 4.1. Photon shot noise simulation ● 5.5 6.0 6.5 7.0 7.5 8.0 8.5 9.0 number of simulated photons (log10 scale) Fig. 4. Correlation between actual and estimated photon numbers. 4.3. Best h estimation (a) (b) (c) Using our image acquisition model, we calibrate the best h for a set of photon density per pixel (dhν ) varying from 20 to 105 . The best h for a fixed dhν is the one that returns the lowest MSE average over 24 images (see Fig.5). Using these results we obtain a relation giving the best h for an estimated dhν (see Fig.6). Using a robust linear, least-square correlation we obtain the following straight line equation: h = 165.2 dhν + 20, 2 again exhibiting excellent correlation (R = 0.9836). 4.4. Image filtering (d) (e) (f) Fig. 3. Shot noise simulations using a sample representative image as a source and a per-pixel Poisson process. From (a) to (c) the mean number of photons per pixel over the whole image are 20, 80, 160. From (d) to (f) the corresponding images with a gain control. Note that there is no simple alternative to this simulation. The resulting Poisson noise is neither additive nor multiplicative, and is not accurately modeled by a simpler distribution (e.g. Gaussian), especially in high-noise, low intensity areas. In Fig.7 we compare the filtering performances of the bilateral filter with those of our proposed adaptive version. For the comparison we simulate sensor images with different dhν , varying from 20 to 320 by doubling at each step the number of simulated photons. We can see that even though the standard bilateral filter has adaptive properties due to its definition it is not adaptive enough to adjust its smoothing behaviour to the noise level. Our proposed noise adaptive bilateral filter has good smoothing behaviour adjustment to noise variations. We can observe that it filters more strongly the image in high noise level while filtering less, favouring edge preservation, in low noise level conditions. Fig. 5. Average MSE calculation over 24 references images Fig. 6. Best h estimation from photon density per pixel 5. Conclusion In this paper, we have proposed a methodology for digital image filtering in the presence of photon shot noise. We have shown that it is important to adapt the level of filtering to the level of noise. In addition, we showed how the bilateral filter can be optimized if the noise level is known. We have proposed a simple and efficient noise level estimation algorithm. This algorithm was shown to exhibit good correlation between estimated and real noise level. By applying this noise level detection method to set the bilateral filter h parameter, we obtain a novel adaptive bilateral filter, which exhibits bestof-class performance at all noise levels. 6. References [1] V.Aurich and J.Weule, “Non-linear gaussian filters performing edge preserving diffusion,” Proceedings of the DAGM Symposium, 1995. [2] Stephen M. Smith and J. Michael Brady, “Susan a new approach to low level image processing,” Int. J. Comput. Vision, vol. 23, no. 1, pp. 45–78, 1997. [3] Carlo Tomasi and Roberto Manduchi, “Bilateral filtering for gray and color images,” pp. 839–846, 1998. (a) (b) (c) (d) (e) (f) (g) (h) (i) (j) (k) (l) Fig. 7. From left to right, the first column shows the noisy images, the second column shows the bilateral filtered images with a h parameter fixed to 100, the third column shows the results obtained with our proposed method. [4] S.Paris and F.Durand, “A fast approximation of the bilateral filter using a signal processing approach,” European conference on computer vision, 2006. [5] Roberto Costantini and Sabine Susstrunk, “Virtual sensor design,” Proc. IS&T/SPIE Electronic Imaging 2004: Sensors and Camera Systems for Scientific, Industrial, and Digital Photography Applications, vol. 5301, pp. 408– 419, 2004. [6] C. Sun, “Fast algorithm for local statistics calculation for n-dimensional images,” Journal of Real-Time Imaging, vol. 7, no. 6, pp. 519–527, December 2001. [7] M. Shinozuka and C.M. Jan, “Digital simulation of random processes and its applications,” Journal of Sound and Vibration, vol. 25-1, pp. 111–128, 1972.