A Bayesian approach for the median filter in image processing

SIGNAL

PROCESSING ELSEVIER

Signal Processing

62 ( 1997) 3033309

A Bayesian approach for the median filter in image processing ’ Giovanni Sebastiani”, Received 31 October

Sebastian0 Stramaglia 1996: revised

16 June 1997

Abstract In this paper we propose the use of the median filter (MF) within the Baycsian framework. This allows us to develop global methods for both image smoothing and image approximation by the MF ‘roots’. A new method for solving the approximation problem is proposed, which is based on stochastic optimization with constraints. Results of the proposed method for both simulated and real binary images are illustrated and compared to results from a known deterministic method. @ 1997 Elsevier Science B.V. Zusammenfassung Es wird der Einsatz des Medianfilters (MF) innerhalb des Bayesschen Ansatzes vorgeschlagen. Hierdurch kiinnen globale Methoden sowohl fur die Bildfilterung als such Rir die Bildapproximation durch die MF-‘Wurzeln’ entwickelt werden. Es wird eine neue Methode zur L&sung des Approximationsproblems vorgeschlagen, die auf der stochastischen Optimierung mit Nebenbedingungen beruht. Ergebnisse der vorgeschlagenen Methode werden sowohl fur simulierte als such echte Binarwertbilder vorgestellt und mit den Ergebnissen einer bekannten deterministischen Methode verglichen. @I 1997 Elsevier Science B.V.

Nous proposons d’utiliser le filtre median (MF) dans un contexte bayesien. Ceci permet de developper des methodes globales a la fois pour I’adoucissement d’images et l’approximation d’images par les ‘racines’ MF. Une methode nouvelle pour la resolution de le probleme de I’approximation est proposee, methode basee sur I’optimisation stochastique avec contraintes. Les resultats de la methode proposee obtenus a la fois pour des images binaires simulees et rcelles sont presentes et compares aux resultats obtenus avec une methode deterministe connue. @ 1997 Elsevier Science B.V. K~ru~~~~d.v:Bayesian inference;

Median filter; Image approximation;

Stochastic optimization

’ This works was made in part at lstituto per le Applicaztoni del Calcolo “M. Picone”, Consiglio Nazionale delle Ricerche, Roma, Italy. * Corresponding author. Address: Istituto per le Applicazioni del Calcolo “M. Picone”, Consiglio Nazionale delle Ricerche. viale del Poltclinico 137. 00161 Roma, Italy; tel.: 396 884 70236; fax: 396 440 4306; e-mail: [email protected]. 0165-16X4/97/$17.00 PII SO 165-16X4(

@ 1997 Elsevier Science B.V. All rights reserved 97)OO 13 I-X

304

G. Srbastiuni,

S. Stramugliu

/ Signal

1. Introduction The median filter (and its derivations) is a useful tool in image processing [ 1,2,5,7, 11, 121. In fact, it can be used, for example, for removing noise from signals while preserving edges. In this context, the median filter (MF) can be used in different ways. In its simplest form, it is a particular type of spatial averaging, that is a very simple smoothing technique which uses the measured values in a small window [ 131. The idea at the basis of this non-linear filter is to provide an estimate of the true grey level at a pixel based on the median value for the observations in a small window (with an odd number of pixels) around that pixel. The median value has the property of being the point whose average difference in absolute value with the observations is minimized. The statistical properties of this estimator are well studied within ‘order statistics’ theory [6, 121. The median estimator works well also when the distribution of the noise has heavy tails. For this reason it is very effective for removing outlier observations. It often works very well in the case of binary noise. However, it performs poorly when the fraction of noisy pixels in the MF window becomes greater than half. Moreover, the performance of the MF is poor when the image is affected by white Gaussian noise [ 131. More complicated schemes are also possible where the MF is used iteratively [3,4, 12, 141. Some of these schemes are based on the property that this filter has ‘fixed points’. These are configurations, usually called ‘roots’, which are not changed by the application of the MF. By consecutively applying the MF to any initial configuration a suitably large number of times, an MF root is obtained. This can be proved theoretically in the one-dimensional case [12, 141. The same property is observed empirically in two dimensions, though without theoretical support. Besides filtering purposes, this property of the MF can be useful for solving approximation problems. This problem consists of approximating an image or a signal by a suitable MF root. This can be relevant for signal compression since MF roots can be efficiently compressed

111. The most natural way of root to any given input image utively the median filter to it procedure usually introduces

associating an MF is to apply consecseveral times. This some unacceptable

Processing

62 (1997)

303-309

distortions. In order to reduce such distortions, a deterministic approach named reduced edge distortion median filtering (REDMF) has been proposed for binary images [4]. Starting from the original input image, local changes of pixel values are repeatedly proposed during a pixel visiting procedure. These changes are accepted if they induced a decrease in the discrepancy between the original image and the MF root obtained by consecutive application of the MF to the actual image. The MF root obtained after the last step of the updating procedure gives an approximation to the input image. Unfortunately, this algorithm suffers from two main drawbacks. First, it contains some logical steps and some parameters which were given for binary images in the case of the MF with a 3 x 3 window, but which have to be determined in other situations. More important, in general this approach will provide a local solution for the problem of minimizing the distortion. When the input image is far from being an MF root, the local solution obtained by the REDMF method can be unsatisfactory. One should then address the problem of finding, among the MF roots, the one which corresponds to the minimum distortion allowed. In this paper, we propose a Bayesian approach (see e.g. [8] and references therein) for using the MF. The image is treated as a Markov random field (MRF). A new type of ‘prior energy’ is proposed for modelling a regularity constraint based on the MF. This approach can be used as a global smoothing technique based on the MF. However, we focus here on the approximation problem. The main goal of using the median filter is then to obtain a good approximation of the original image by an MF root. Small differences between the original and the processed images are allowed, without losing important features. An approximated version of the simulated annealing (SA) algorithm [9] is used here for solving globally this constrained optimization problem. Should the result of the approximated SA procedure happen to be different from an MF root, the subsequent application of the REDMF method to it will provide an MF root. The proposed method can be applied without changes to different situations like the MF based on different windows. As already observed, median filtering performs well in the presence of binary noise (which is the only type possible in binary images), while it performs poorly when the noise is not binary,

e.g. Gaussian (which is most common in grey level images) [ 131. We, therefore, concentrated our study on binary images.

2. Description of the method Let y be an original N x N digital image. We introduce an image x which can be either the MF root approximating y or the non-degraded version of y (when using the approach for smoothing). We model the image x as an MRF with probability function n(x) given by the following Gibbs distribution:

where H(x) = Ilx-9(x))12, the norm is the quadratic one, and the summation is over all possible digital values in an image. The local operator 9 is the MF spatial averaging, which transforms the image x by substituting the value at any pixel by the median value in the MF window centered at the pixel. The above choice for n(x) gives the largest probability to the MF roots since H(x) is equal to zero only in correspondence to these configurations. In order to discriminate between two different MF roots, we add to H(x) a ‘fidelity’ term L(x,y) = I/Y - x/I? Pena 1’izm g configurations x which are not ‘close’ to the original image y. The considered function f?(x) will then be fi(x)=L(x,y) + X(x), where i. is a positive parameter. The above written functional can be introduced into Eq. (I ) to compute probabilities. For a given value of i., the mean value or the mode of the probability function n(x) can be considered. Approximations of these quantities can be obtained by sampling from n(x) via Markov chain Monte Carlo (MCMC) algorithms [lo]. This will provide configurations that approach MF roots as i is chosen larger and larger. On the other hand, when 3. is close to zero, the configurations approach the original image y. With a suitable choice for & the obtained configuration could be used for smoothing purposes. Recall that here we focus on the approximation problem which consists of finding the MF root which is closest toy. We will now see how the stochastic approach described above can be used to solve this problem. We observe that the constraint for x to be an MF root is given mathematically by the relation H(x) = 0.

The approximation problem consists of finding the x which minimizes the fidelity functional L(x, y) among those verifying the constraint H(x) = 0. This problem can be solved in general for any positive bounded functionals L and H by using stochastic optimization with constraints [8]. In fact, let n,(x) be the probability measure obtained by introducing into Eq. (1) the functional H,(x)=fi(t)(L(x,y) + i.(t)H(x)), where P(t)>O. For any fixed t, we can consider a suitable irreducible and aperiodic Markov chain which has n,(x) as its unique invariant measure [8]. Starting from an arbitrary initial point x0, the chain is generated by single pixel updating within a visiting procedure (random, raster scan). The new grey level value X, at pixel i is randomly drawn from the conditional probability of x, given the values of x in all other pixels. Due to the MRF nature of x, this conditional probability depends only on the values of x at pixel i and at pixels of a small neighbourhood of pixel i. This leads to a small computation burden for pixel updating. It can be proved that the limit of this chain follows the distribution n,(x) independently of the distribution from which the initial point of the chain is drawn [9]. Instead, let us make /3(t), A(t) + x under the condition lim ,_&l(t)i(t)/ log t = c < 1/,Y2/&, where H,,,= max, H(x). Then, the chain will converge to one of the configurations x satisfying the constraint H(x) = 0, and for which the fidelity functional L(x, y) has its absolute minimum. This solves the approximation problem. The conditions stated above, which guarantee convergence cannot be fulfilled during implementation because of computing limitations. In practice, we adopt an approximated schedule with p(t) and i(t) growing linearly, as described in Section 3. Although there is no theoretical guarantee that the approximated algorithm will converge to an MF root, this usually happens in practice. However, when some changes appear after applying the MF to the result of the proposed method, they are usually few. In that case, an MF root is obtained by applying the REDMF method to the result of the proposed method. This provides a very good approximation of the input image by an MF root. In fact, the result of the proposed method usually has a very small discrepancy from the input image. Moreover, it is very close to being an MF root, so that the REDMF method works very well.

306

G. Sebastiani,

S. Stramagliu / Signal Processing

3. Results and discussion In this section we describe and discuss some results obtained from the proposed stochastic method when applied to artificial or real binary images. The results are also compared to those obtained by the deterministic REDMF method [4]. Since the REDMF algorithm has been described only for the MF with 3 x 3 window (second order neighbours), we first implemented our algorithm with a 3 x 3 window. For each experiment, over 10,000 image ‘sweeps’ (single updates of all image pixels) were used for the proposed method. Pixels were visited following a raster scan procedure. The initial configuration was chosen randomly. The parameters p and 3, were increased linearly in the sweep index in the intervals [l, 71 and [0, lo], respectively. The domains for /3 and i, were chosen by trial and error for the error-free synthetic image presented here. We then used those values without adjustments to process many other images of different objects and types, and also based on different MF windows. We note that the adoption in the proposed method of a linear variation of the parameters fl and L was imposed only by computation time limitations. Instead, the theoretical variation would have guaranteed algorithm convergence. This situation is completely different from that of the REDMF method, where a list of the values of the existing parameters is simply given, without any alternative theoretical rule to choose them, and only for the 3 x 3 window case. The regular recursive scheme was adopted in the REDMF method to filter the image to an MF root at each step that this was required. We implemented the two methods described above by using the Fortran 90 language on a serial computer Alpha 2100 (Digital) with a clock frequency at 275 MHz. Fig. l(a) shows an artificial 64 x 64 image. The MF root in Fig. l(b) is the result obtained from the proposed method. As can be seen, no significant part of the original image was lost after the processing. This was not the case for the root obtained from the REDMF method, as shown in Fig. l(c). We note that in areas where pixel values in the original image are unchanged by the MF as spatial averaging, both methods obviously induce no distortions. Small distortions from both methods are induced when ‘isolated’ (on scale of the MF window) pixels are present, whose

62 (1997) 303-309

value is changed by the MF. On the other hand, when such pixels are present sufficiently close to each other, the local procedure for pixel updating of the REDMF method does not always allow us to recover the MF root closest to the original image. Distortions in the REDMF image can then appear in aggregated form, so that macroscopic structures are altered. In fact, in this example the three small squares and a ‘bridge’ in the original image disappeared in the REDMF image. This did not happen for the proposed method as seen in Fig. l(b). In fact, the proposed method replaced, for example, the small squares with the crosses, that correspond to the MF roots that have minimum distortion from the squares. This different behaviour is essentially due to the fact that we are using a global optimization procedure, while the REDMF method is working locally. The difference (measured by the fraction of misclassified pixels) between the original image and the MF roots were 0.98% and 2.0% for the proposed and the REDMF methods, respectively. It appears that the proposed method performed better than the REDMF method both visually and quantitatively. The stability of the proposed method with respect to presence of noise in the original image was also assessed. We tested the method on an image obtained by adding binary symmetric channel noise with probability 0.02 to the one in Fig. l(a), as shown in Fig. l(d). By applying the proposed method to the noisy image, we obtained the MF root shown in Fig. l(e). The percentage distortion of this MF root from the original image was 1.3%. This indicates that the method is robust to the presence of noise, which only slightly degraded the processed image. We also applied the method to real images. As an example, Fig. 2(a) shows a 54 x 273 hand written binary image of the words ‘median filter’. The binary image was obtained by scanner digitalization and successive thresholding. Moreover, in this case, the proposed method produced an MF root (Fig. 2(b)) in good agreement with the original one. On the other hand, the REDMF method performed very poorly, as seen in Fig. 2(c). In another example, a 200 x 360 binary image was produced as described above starting from a photograph of rocks obtained by electron microscopy (Fig. 3(a)). The results from the proposed and the REDMF methods are presented in Figs. 3(b) and (c),

G. Sebastiani,

S. Swamaglia

/ Signal

Processing

307

Fig. I. An artificial image. (a) The original image; (b) the MF root obtained by the proposed method; (c) the MF root obtained by the REDMF method; (d) the noisy version of the image in (a); (e) the MF root obtained by applying the proposed method to the image in (d).

Fig. 2. A real image representing the handwritten words ‘median filter’. (a) The original image; (b) the MF root obtained by the proposed method; (c) the MF root obtained by the REDMF method.

respectively. We note that the REDMF method often fails in separating two white objects which are close to each other in this image while merging them into a single region. This did not occur for the proposed method, except in one pair of objects. We observe that for some applications a relevant quantity to be estimated

is the distribution of the size of the particles (objects) in the picture. It follows that in such applications the possibility that distinct objects are merged after the processing should be minimized. We also present in Fig. 3(d) our algorithm implemented with a ‘cross’ window (first-order neighbours). We observe that more details are preserved when using a smaller window. This fact is expected for a method that provides an MF root not far from the input image. In fact, space resolution in the best approximation of an image by an MF root increases as the MF window size decreases. This is due to the following two facts [ 121. First, the size of the so called ‘minroot’, that is the minimal structure that is invariant when applying the MF to it, decreases as the window size decreases. For example, for binary images the minroot of the cross MF window is the 2 x 2 window, while the minroot of the 3 x 3 MF window is the cross with arms two pixels wide. Secondly, structures whose sizes are smaller than the minroot size are destroyed by applying the MF and therefore cannot appear in the MF root approximating the input image. We could not compare this last result to the corresponding one from the REDMF method with the cross window, since this required knowledge of some

308

G. Sebastiani,

S. Stramaglia / Signal Processing

62 (1997) 303-309

logical steps and of some parameters, which were given by the authors only for the 3 x 3 MF window case. Finally, we note that the computation time for the proposed method was larger than for the REDMF method when the original image was not too far from being an MF root. This happened, for example, in the artificial image presented here. In this case, the ‘CPU times for the proposed and the REDMF methods were about 200 and 66 seconds, respectively. However, when this condition was not fulfilled, the computation time for the proposed method was comparable if not smaller than that for the REDMF method. For example, in the case of the ‘median filter’ image included in this paper, the ‘CPU’ times for the proposed and the REDMF methods were about 770 and 1270 seconds, respectively.

4. Conclusions (b)

The use of the MF within the Bayesian framework was proposed. A new stochastic method for image approximation by MF roots was developed. The proposed method can be used without changes in different situations, like MF windows of any kind. The application of the method was successful for both artificial and real binary images. The method performed better than a related deterministic method described in the literature.

Acknowledgements The authors are thankful to F. Godtliebsen for his suggestions which improved our presentation and to G. Dellino for providing the electronic microscopy data.

References (4 Fig. 3. A real image representing rocks. (a) The original image; (b) the MF root obtained by the proposed method; (c) the MF root obtained by the REDMF method; (d) the MF root obtained by the proposed method implemented with a cross window.

[I] G.R. Arce, N.C. Gallagher, BTC image coding using median filter roots, IEEE Trans. Commun. COM-31 (June 1983) 784-793. [2] D.G. Bailey, R.M. Hodgson, S.J. McNeil], Local filters in digital image processing, Proc. 21 st IEEE National Electronic Conf., Christchurch, New Zealand, 1984, pp. 95-100

[3] W.W.

Boles,

dimensional extraction. 1988) [4] W.W.

Kanefsky, filtering

IEEE

M.

Simaan,

algorithms

Trans. Circuits

Recursive

two-

Boles,

Systems CAS-35

(October

21

[5] R.T. Chin.

Kanefsky,

Yeh,

noise

A

reduced

for binary images. Signal

Quantitative

smoothing

techniques.

( 1983) Wiley,

[7] B.R. Frleden. A new restoring

of some edge

Computer

Vision

67-9 I. 1970.

( 1976)

280.-283. [8] D. Geman, P.L.

Hennequin.

de Salni-Flour Vol.

tributions.

(Ed.).

XVII-1988.

1427. Springer,

[9] S. Geman.

]

R.M.

Hodgson,

rank filters.

D.

Berlin.

Geman,

and the

&ole

D’!%

Lecture

Notes

ProbabilitO

in Mathematics.

1990. pp. 113-193.

Stochastic

Bayesian

de

in imaging.

relaxation,

Huang

restoration

of

images,

D.G.

dihIEEE

method< using Markov

Biometrlka

Bailey,

Properties.

M.J.

( I ) ( 1970)

57

Naylor.

implementations

Image Vision (Ed.).

Comput.

I

Two-DimenGonal and Median

Phy-sits. Vol. 42. Springer, 131 .A.K.

Jain.

A.L.M.

Ng.

and applications

(I )

(February

ot

1985)

Fundamentals Englewood

Berlin. of

Digital

Signal

I98 I, Chapters Digital

ClifTs.

Pro-

Filters. Topics in Applied I. 5 and 6.

Image

Processing,

1980.

Chapter

NJ.

7.

pp. 246-249. 141 P.D.

Wendt.

vergencc

E.J. Coyle,

properties

Syst. CAS-33 Gibbs

1984)

3-14.

Prentice-Hall. Random fields and inverse problems

applications.

cessing II: Transforms

for the preferential

of edge gradients. J. Opt. Sot. Amer. 66

6 (6) (November

97-109.

121 T.S.

New York.

algorithm

and their

S.J. McNeill.

evaluation

Intell.

IO] W.K. Hastings, Monte Carlo sampling

edge

1990) 37-47.

Order Statistics,

enhancement

Simaan,

algorithm

Image Process. 23

David.

M.

( I ) (September

C.L.

preserving Graphics

M.

Machine

721-741. chains

median filtering

Processing

in:

Trans. Pattern Anal.

for fast image root

1323-1326.

distortion

[6] HA.

M. median

(March

N.C.

of median

Gallagher, filters,

1986) 276-286.

IEEE

Jr.. Some

con-

Trans. Circuits