22 Order statistics in image processing

N. Balakrishnan and C. R. Rao, eds., Handbookof Statistics, Vol. 17 © 1998ElsevierScienceB.V. All rights reserved.

Z,L

Order Statistics in Image Processing

S c o t t T. A c t o n a n d A l a n C. B o v i k

I. Introduction

In the age of the home computer and desktop workstation, the processing, analysis, storage and display of visual information as digital imagery is becoming commonplace. High definition television (HDTV), direct TV, photo CDs, the World-Wide Web and video teleconferencing are all examples of emerging technologies that utilize digital images. This technological explosion has also been accompanied by many advances in the processing of digital media. Thus, digital image processing has become one of the fastest-growing areas of research and development in the electrical and computational sciences. Current areas of significant research in image processing include the development of theory and algorithms for the acquisition, display, decomposition, reconstruction, compression, enhancement, restoration and analysis of digital image data. The environments in which digital image data are captured are often subject to unpredictable noise levels, interference and modeling uncertainties. Thus, statistical techniques in general have been widely adopted in digital image processing and numerous papers have elaborated optimization criteria for performing image processing tasks. Not surprisingly, in practical applications, the lack of knowledge of signal contamination suggests the use of robust techniques to accomplish many of these tasks, and so, order statistics have found an important place in the image processing literature. Although the statistical properties of local order statistic-based estimators are well-understood, and their application to signal and image processing might be guessed at, there are many serendipitous non-statistical attributes of order statistic methods that make them particularly appealing. Some of these attributes are so striking that order statistic based methods have come to dominate certain areas of computerized image processing. It is hoped that making these applications accessible to statisticians, applied mathematicians, non-electrical engineers, etc., will increase the cross-fertilization of new concepts. Thus, this chapter seeks to summarize several of the most significant applications of order statistics in the area of digital image processing. It is appropriate at this point to give some background into the problem of digital image processing, and the overall motivation for using order statistics in

603

604

S. T. Acton and A. C. Bovik

image processing problems. Digital images may be obtained by sampling and quantizing a continuous image, whether color or black-and-white. The continuous image may be represented as a continuous electrical waveform, such as a television signal, or as a physical negative image, positive slide, or photograph (strictly speaking, a photograph is quantized at the level of the grains in the emulsion; usually, however, this is much finer than the sampling resolution). Digital images may also be acquired directly, as from a digital camera, for example; such instruments are now becoming widespread and extremely inexpensive. Often, the image obtained will be corrupted by noise or cluttered with unessential detail. The noise may arise from electrical noise, from imperfections in the imaging environment (ranging from faulty sensors to dust on the camera lens!), or from channel errors arising in the transmission of the digital data. The techniques of digital signal processing (DSP) offer powerful methods to improve and analyze these data, and have found application in the area of image processing for over three decades. Most basic among these tools is the digitalfilter, which most generically, is a transformation of a digital signal or image into another signal or image, according to some prescription. However, most of the existing theory, and certainly, practical implementations have involved the theory and design of linear digitalfilters. Linear filters are, of course, transformations that obey the linear superposition property, examples of which can be found in many aspects of everyday life. Automobiles are equipped with mufflers (acoustic filters) and shock absorbers (mechanical filters); hi-fi systems are supplied with electrical filters (now often of the digital variety) that provide a clear signal output. Such filters can be implemented (or modeled) by linear convolution, and can be described by a frequency response. Thus, linear filtering may always be viewed as a process which amplifies, attenuates, and phase shifts the individual frequency components of the input signal. The design of linear filters thus always implies some specification based on frequency selection. Linear filters can be expeditiously implemented in the frequency domain via the Cooley-Tukey Fast Fourier Transform (FFT) techniques, or other recent fast algorithms. Indeed, the initial euphoria surrounding signal processing in the 1960s may be attributed to the efficacy and the straightforward implementation of FFT-based linear filters. Despite their many advantages, linear filtering approaches suffer from certain limitations which have emerged in certain application domains, and which prove to be difficult to surmount without abandoning the linear approach entirely. For example, a common task in image processing is the removal of broadband additive noise (e.g., white noise), which occurs in many image signals, such as radar data, electron micrographs, tomographic images, etc. When the noise is broadband, it is inevitable that the signal and noise frequencies will overlap, often in regions of the spectrum that support key pieces of signal information. Linear filtering methods have no ability to distinguish between noise and signal components occupying the same parts of the spectrum, and so, the process of noise reduction can also result in the distortion of important signal properties such as image edges and image detail.

Order statistics in image processing

605

Thus, nonlinear image filtering, and as we will see nonlinear image processing paradigms in general, have become an area of energetic research. Numerous techniques have been proposed, implemented, demonstrated, and argued, but in only one area of nonlinear image processing has there emerged a substantive discipline based on mathematical rigor, general applicability, and repeatable performance. That area may be accurately described as order statistic-based image processing, although the wealth of research that has emerged has revealed a wide variety of new nonlinear methods that are not intrinsically involved with order statistics, but which are, nevertheless, direct descendants of order statisticbased methods. In this overview, we will highlight several established developments in the area of order statistic-based image processing. In addition, a number of recent innovations and extensions to the basic theory will be illumined. Throughout, the practical application of these techniques to digital signal processing, and image processing in particular, will be emphasized. At this point in time, a number of image processing devices, such as the median filter, the order statistic filters, and many related techniques are no longer novelty items, but instead, are standard image processing tools. The goal of this contribution is to give an overview of the role that order statistics have played in this. The plan is not to be exhaustive, but instead, to give suggestive results while maintaining a modicum of analysis. Because of the disparity between the notations usually used by statisticians and by image processing engineers, we shall adopt a simple format and exposition in an effort to maximize communication.

2. Order statistic filters

An M-dimensional digital image will be denoted {Xn : n C V C ZM}, where Z is the set of integers, and where the domain V is discrete and countable (and in practice, finite): typically an M-dimensional integer interval. Such digital images are also quantized to have discrete, finite ranges; computer processing requires this, although the analysis does not always make this explicit. The elements xn represent samples of the image, called pixels (picture elements). Usually, this equates to samples of optical image intensity, or brightness (which is positivevalued), but it may also represent other attributes of sensed radiation. Indeed, the diversity of image types is as diverse as the different bands of sensed radiation, which may be electromagnetic, sonic, or atomic, and also dependent on the types of sensors used to capture the data. Some images may even be complex-valued, as in the case of synthetic aperture radar. However, for simplicity of exposition, we will assume that the samples Xn are real-valued samples of optical images, and all of the image processing examples given will be standard optical images; however, the techniques summarized in this chapter can be, and have been, applied to an amazing diversity of image types. Two-dimensional digital images (the most common variety, representing samples of static, spatial imagery) are usually represented as a two-dimensional

606

s. T. Acton and A. C. Bovik

matrix of values, {xi,j}, where i indexes the image row and j the image column. However, in order to maintain generality, we will utilize the M-dimensional form x,. There are two reasons for this. First, many image signals are of a higher dimensionality, the most obvious example being digital video, where there are typically two spatial indices and one time index. There are also images naturally having three spatial dimensions, such as tomographic data or optical sectioning microscopic images. The second reason we will adopt the general M-dimensional form is that most order statistic-based processing paradigms do not make use of the spatial ordering of the data beyond gathering the data locally; instead, the algebraic ordering is considered (prior to weighting, processing, etc.). Regardless of dimensionality, a digital filter transforms a signal {x, : n E V C_ Z M} (hereafter simply {xn}) into another, filtered signal {Yn}, usually of the same dimensionality and domain, with the intention that the result will be somehow improved: less noisy, for example, or with certain features enhanced or emphasized. Median and rank-order filters

The median filter (Tukey 1971; Tukey 1974) is a simple nonlinear device which does precisely what linear filtering paradigms could not accomplish - remove high-frequency noise from signals without disturbing certain high-frequency informational components - specifically, sharp, sustained transitions and trends, such as edges at the boundaries of objects, or other details representing minute, yet significant structure in the image signal. And more - it had the capability to eradicate impulsive noise occurrences with amazing efficiency. And best of all, it was conceptually trivial and easy to implement - all one had to do was define a moving window, and replace incoming signal points with the computed sample medians. Shortly after the conception of the median filter, a number of DSP engineers immediately foresaw that it might prove to be useful for such significant application tasks as speech signal filtering (Rabiner et al., 1975) and image filtering (Frieden, 1976). However, at that time there was no DSP theory for the median filter, and while there did exist a vast statistical literature on the sample median and other order statistics, it was both written in a different technical language, and dealt mostly with static population samples, and very little with signals. Since then, a rich DSP theory of order statistics has begun to emerge, and which may be regarded as the DSP engineer's contributions to the field of order statistics. A significant part of this contribution has been made by engineers concerned with image signals, which are strongly characterized by the type of information that the median filter tends to favor. Given an M-dimensional image {x.}, the median filter requires the definition of an M-dimensional moving window with a fixed number N of elements, to compute each individual output value. A window is really a pre-specified geometric or temporal law for collecting image samples around each coordinate, in order to perform the median or other operation on them. Once the samples are collected within the window, their relative positions in space (within the window) are discarded, and so it suffices to define a one-dimensional vector to contain

607

Ol'der statistics in image processing

them. Therefore, at image coordinate n, define the windowed set of pixel values by the one-dimensional vector x. =

(1)

This notation indicates that there are N image pixels, indexed 1 , . . . , N , that have been captured from the vicinity of image coordinate n. It will always be assumed that the window covers an odd number N = 2m + 1 of image samples. This is easily justified: at each image coordinate n the filter window is centered directly over the sample to be created at the same coordinate in the output image; the window is usually desired to be (even) symmetric along each dimension; and the window ordinarily contains the current pixel over which the window is centered. Median filter windows

The windowing law used to collect spatial image samples operates according to a geometric structure such as a square-shaped, cross-shaped, hexagonally-shaped, or other such window, several of which are depicted in Figure 1. Usually, windows having a nearly circular geometry are favored, since this encourages isotropic (direction-insensitive) processing of the data (the ROW and C O L U M N filter geometries being an exception to be explained shortly). Of course, it is sometimes desirable to encourage directional processing, in which case a more

m D 3x3 SQUARE

m m 5x5 SQUARE

ikdl

m 3x3 CROSS

/

m 5x5 CROSS

um

ml 3x3 XSHAPE

5x5 XSHAPE

lx3 ROW and 3xl COLUMN

lx5 ROW and 5xl COLUMN

Fig. 1. Typical window geometries for the median filter and other order statistic filters. The white square denotes the center pixel covered by the window.

608

S. T. Acton and A. C. Bovik

oblate window may be selected. The windowing law is defined by a set of vector s h i f t s B = { b i , . . . , bN}, (usually centered around the zero vector shift 0) that determine the coordinates of the image samples that will be captured into the windowed set. With respect to a given sample x,, the window B may be considered an operator (a kernel) such that (Xn__bl,... ,Xrl_cbN)T = XtT. Thus, in Fig. 1, the window B = ( ( - 1 , - 1 ) , ( - 1 , 0 ) , ( - 1 , 1), ( 0 , - 1 ) , (0,0), (0, 1),(1,-1), (1,0),(1, 1)), and at each coordinate n = ( i , j ) the windowed set associated with B is coordinate

Xn = (Xi- 1, j - 1, Xi- 1,j ~X i - I ,j+ 1 : Xi,j- 1, Xi~j, Xi,j+ 1, Xi+ 1, j - 1 ~Xi+ 1d' Xi+ 1,j+ 1) T .

Before actually defining the simpler filters, one other detail needs explanation, When the index n of window of a nonlinear filter is brought next to or near the edge of the image, part of the window will overlap with "empty space"; viz., certain elements of x. will be undefined. This is depicted in Figure 2 for a 3 x 3 S Q U A R E filter. Although there is not strong theoretical justification for it (and other methods are used), the usual convention that is followed in defining the windowing operation near the borders of the image is to (conceptually) fill the "empty" window slots with the value of the nearest well-defined image pixel nearest. This is called p i x e l r e p l i c a t i o n . This is usually regarded as more satisfactory than filling with zero values, since it implies that the scene being imaged falls to a zero intensity beyond the border, which is generally not the case. Another approach is to r e f l e c t pixel values, i.e., define the windowed pixels beyond

0

•

0

0

Fig. 2.Depiction of a 3 x 3 SQUARE window overlapping the borders of an image at different locations.

Order statistics in image processing

609

the image border by symmetrically reflecting pixels across the axis defined by the image border. This often has the advantage of preserving the continuity of patterns, but it is also harder to program.

Median filter definition and examples With this notation in hand, it becomes possible to define and develop a large class of nonlinear filters. Given a windowed set x,, of image samples collected at and near coordinate n, define the vector of order statistics of x, by x(.) = {x(1):.,...,x(N):.} r ,

(2)

where x(1):. =_< • • • _< X(N):,. The output of the median filter is then easily defined: let {y.} = med{x.} when Yn = X ( m + l ) : n

•

(3)

As an example of the median filter's ability to eliminate impulse noise, see Fig. 3. Figure 3a displays a digitized optical image of a tower on the University of Texas at Austin campus. The image has been digitized to a spatial resolution of 256 x 256 pixels (65,536 pixels total), with each pixel's intensity quantized to one of 256 possible gray levels (8 bits), indexed from 0 (black) to 255 (white). The image in Fig. 3(b) contains "salt and pepper" impulse noise - 25% of the image brightness values have been randomly changed to extreme white or black values (with equal probability). Figures 3(c), 3(d), and 3(e) show the result of median filtering the noisy image in Fig. 3(b) using (square) windows of sizes 3 x 3, 5 × 5, and 7 x 7 (covering 9, 25, and 49 pixels), respectively. The quality of the enhanced images in Figs 3(c)-3(e) is quite variable - there is a strong dependence in general between median filter results and filter geometry, as might be expected. The 3 x 3 filter was not able to remove clustered outliers very effectively - although the overall density of impulses was reduced, others groupings of outlying values were merged into larger structures, which could lead to misinterpretation of their nature. By contrast, using a larger 5 x 5 filter window provided superb noise suppression - certainly better than attainable with any linear-type filter, while not degrading the informational part of the image much - the edges, or object boundaries, and much of the detail remains sharp. This result can definitely be interpreted as an enhancement of the image. Taking a yet larger filter window, however, as evidenced by the result attained using the 7 x 7 window, can lead to excessive destruction of some of the useful information. There is clearly a tradeoff between noise suppression and information retention in the selection of an effective median filter window for image smoothing.

Median filter roots Prior to 1980, little work was done on characterizing the temporal or spatial properties of the median filter - as with all nonlinear filters, difficulties in analysis often lead to intractable problems in time- or space signal processing. Indeed, the only statistical treatment on the subject prior to 1980 that we have been able to

610

S. T. Acton and A. C. Bovik

(a)

(b)

(c)

(d)

(e)

(13

Order statistics in image processing

611

find that actually incorporates time- or space-dependent information is H. A. David's work characterizing the correlations of ranges from overlapping samples (David 1955). This work has since been extended towards developing the spectral properties (second-order distribution and selected moment properties) of stationary signals that have been median filtered or filtered by other order statistic-related filters (Bovik and Restrepo 1987; Kuhlman and Wise 1981). The first significant characterizations of the non-static properties of the median filter demonstrated that certain interesting one-dimensional signals, called root signals (or fixed points) are invariant to one-dimensional median filtering (Gallagher and Wise 1981; Tyan 1981; Longbotham and Bovik 1989). While we do not plan to spend much time in discussion of one-dimensional signals (they are discussed at great length elsewhere in this volume), some mention of the basic root signal properties will enhance the discussion of two-dimensional signals. The root signals are characterized by a novel measure of signal smoothness, termed local monotonicity. A one-dimensional signal {x, } is locally monotonic of degree d or LOMO-(d), if every subsequence of {x,} of d successive elements forms a monotonic (increasing, decreasing, or constant) sequence. THEOREM 1. Tyan (1981), Longbotham and dimensional signal {xn} contains at least one length m + 1. Then the output of a y, = X(m+l):~ = Xm+~:, = X~ for every n if and

Bovik (1989): Suppose that the onemonotonic segment (xk,... ,x~+m) of length N = 2 m + 1 median filter only if {x~} is LOMO-(m + 2).

This unexpected result gives deeper insights into the type of signal that the median filter "prefers," and it also suggests the possibility of a limited eigenfunction-like analysis for median filter design and performance study. Even more surprising is the fact that repeated median filtering reduces any signal to a root: repeated application of the median to most finite-length signals (only excluding oscillating bi-valued signals) results in convergence to a root signal in a finite number of iterations: THEOREM 2. Gallagher and Wise (1981), Tyan (1981): A one-dimensional median filter with window size N = 2m + 1 will reduce a length-L signal {xn} to a root signal that is LOMO-(m + 2) in at most (L - 2)/2 repeated passes. Some mention needs to be made at this point about signal duration. In Theorem 1, it is assumed that the signal may be doubly infinite in duration, although it have a finite support. In Theorem 2, the signal is assumed finite in extent (and undefined beyond) - and so there must be some mechanism for defining the filter operation

Fig. 3. (a) Original "Tower" image. (b) "Tower" corrupted by 25% salt-and-pepper noise. (c) 3 x 3 SQUARE median filter applied to noisy image. (d) 5 x 5 SQUARE median filter applied to noisy image. (e) 7 x 7 SQUARE median filter applied to noisy image. (f) 50 iterations of 3 x 3 SQUARE median filter applied to noisy image

612

S. T. Acton and A. C. Bovik

(so that the window is filled with samples near the signal endpoints). The agreed upon method of doing this is to "pad" the signal, by appending m samples at the beginning of the signal (all equal to the first sample), and m samples to the end of the signal (all equal to the last sample). This has the added benefit of introducing the requisite monotonic segment of length m + 1 of Theorem 1. In fact, for finitelength signals with padding, Theorem 1 can be modified. In the context of image filtering, this process becomes identical to pixel replication. Both of these results hold only for signals of one dimension; similar results for true two-dimensional windows, such as SQUARE, CROSS and X-SHAPE in Fig. l, have not been developed. Nevertheless, understanding the behavior of the one-dimensional median filter, and its root structures, does yield insights into its efficacy in image processing applications. In image processing parlance, the median filter is frequently referred to as an "edge-preserving, impulse-suppressing" image smoothing device. For one-dimensional signals, the preceding Theorems make clear that the median filter does not disturb signal trends - even sharp ones, and of course, the sample median possesses high efficiency for estimating data contaminated by outliers or by heavy-tailed noise influences. We will not spend much time on the statistical characterization of the median filter here, other than reviewing some distribution formulae for two-dimensional windows. The reason for this is that, although statistical considerations are equally important for image processing applications as to standard statistical applications, most of them will elsewhere be available in this volume, or already familiar to the reader. As mentioned above, there has been no theory developed that describes the root-signal structures of median filters defined with two-dimensional windows, such as those depicted in Fig. I. Indeed, it is not at all clear that a well-defined class of such root signals exist, although it is trivial to construct examples of twodimensional root signals for any median filter window (for example, a two-dimensional signal that is constant along one dimensional, and locally monotonic of sufficient order along the direction orthogonal). Part of the reason for this is that there is no clear definition for two-dimensional local monotonicity, since it is the orientation of each image path along which the pixel values are to be locally monotonic must be considered.

Separable median filter Actually, there is a type of two-dimensional median filter that does admit a root signal analysis, although the filter windows that define it are not truly two-dimensional. The so-called separable median Jilter is defined by filtering the image row-by-row with a horizontally-oriented filter window, then filtering the result column-by-column with a vertically-oriented filter. Figure 1 depicts length-three and length-five ROW and C O L U M N windows that define the separable median filter. The term "separable" is not accurate, since the so-called separable filter is not equivalent to a true two-dimensional filter result (i.e., not equivalent to the CROSS filter). Also, row filtering followed by column filtering will not generally yield the exact same result as column filtering followed by row filtering, although the two results will generally be qualitatively nearly identical.

Order statistics in image processing

613

The roots of the separable median filter are described by Nodes and Gallagher (1983), who found that the two-dimensional roots are images that must be roots of both the row and column filters individually (locally monotonic along both directions), except for certain instances where there may be an image region over which the image is strictly bi-valued, and oscillates in a certain way in that region (Nodes and Gallagher 1983). However, such binary oscillating regions are so rare in digital imagery that the bi-directional root property practically characterizes the two-dimensional roots of the separable median filter. Stated another way, the practical roots of a separable median that utilizes windows of span 1 x (2m + 1) and (2m + 1) x 1 are images whose rows and columns are all LOMO-(m + 2). This property approximates many naturally-occurring images, and often is taken as a simple and natural definition of two-dimensional local monotonicity. Nodes and Gallagher also are able to show that repeated application of a separable median filter to an image of finite duration will nearly always eventually reduce the image to a root in a finite number of applications. The exceptions are again, certain bi-valued oscillating images that are not encountered in practice. Further root image analysis for recursive modifications of the separable median filter have been supplied by McLoughlin and Arce (1987). Two-Dimensional roots

However, there has as yet been no generalized root image theory developed for true (non-separable) two- and higher-dimensional signals, nor is there evidence that repeated two-dimensional median filter applications must converge to a root signal. Thus, M-dimensional root signal theories remain the subject of intense inquiry. Not much progress has been reported, although for strictly bi-valued images, Heijmans (1994a) has introduced a modified median filter that will converge to a fixed signal. Certainly, part of the problem lies in the fact that the broad diversity of possible two-dimensional windows makes simple root-image characterization unlikely. Further, experiments do not show that two-dimensional median filters yield root images upon repeated application. For example, Fig. 3(f) depicts the divergent behavior of the median filter for digital imagery. After 50 iterations of the two-dimensional median filter, blurring of the image continues to evolve, the image modifying slightly at least with each application. Regardless of the existence or otherwise of two-dimensional root images, the median filter demonstrably smoothes image information quite effectively is such a way that information-bearing sharp transitions, or "edges", are maintained. Image edges contain a wealth of information about the objects contained in the image: their locations, shapes, sizes, and contrast relative to the background. Much of image analysis theory and algorithm development relies upon first detecting the edges in a digital image. By contrast, low-pass linear filtering invariably blurs edges and promotes inter-region smoothing, thus reducing much of the key information available in the edge data.

S. T. Acton and A. C. Bovik

614

Streaking and blotching in median filtered images Because of its amicable properties for image processing, the median filter has enjoyed tremendously widespread use in an amazing diversity of image processing ~ystems, including, for example, commercial biomedical image processing systems, medical radiographic systems, military imaging hardware, etc. However, the median filter is not a panacea! A significant drawback associated with median filtering are the phenomena known as streaking and blotching (Bovik 1987). When a median filter is iteratively applied to an image, or when a filter is used that has a large window span, significant runs of constant or nearly constant values can create streaking artifacts in the output signal. These are a consequence of the median filter's tendency to create root-like signals (see Theorem 2). Bovik (1987) supplies an extensive analysis of the streaking and blotching effects of the median filter. The statistics of median filter streak lengths for one-dimensional signals can be arrived at straightforwardly, although the analysis of two-dimensional streaks (blotches) is rather involved. For the case of one-dimensional filters, suppose that the one-dimensional signal xn is a sequence of independent and identically distributed (i.i.d.) random variables with common arbitrary continuous probability distribution function F. As before, denote the output of a length N = 2m + 1 median filter by yn = X(m+l):,; define a streak of length L in y~to be an event of the form

{Yn-1 7g Yn = Yn+l . . . . .

Yn+L-I • Y,+L} ,

(4)

that is, a constant-valued output subsequence of duration exactly L. Given the initial coordinate n, L is a discrete random variable with a finite range space (since, under the assumption of continuity of F, nonoverlapping windows will yield equal outputs with zero probability). The probability mass function of the streak length is given by 2m÷ 1

PL(fl~)=Pr{Z=2}-- m+-l-[Qm(fl~+2)+Qm(.~)-2Qm(fl~+l)]

(5)

for 1 < 2 < 2m + 1 and PL(2) = 0 elsewhere, where, defining rl = m a x ( l , 2 - m), r2 = min(2, m + 1)

Qm(;O=

er{yn

=

.....

2m-2+2~

2m + 2

Y.+L-1 }

(2r-I)km-r+lJ r=rl

(6)

{2m+a-?~

k m+r 1 J

for 1 < 2 < 2m + 1 and Qm(2) = 0 elsewhere. In the above expression, there is no dependence on the window position n, since the input {x, }, and hence, the output {y,}, is stationary. Note also that there is no dependence on the input distribution F, other than the property of continuity. Thus, median filter streaks tend to occur with lengths independent of the input distribution. It can be shown that for large windows spans N = 2m + 1, PL(2) --+ 2 -~. The expected streak length is

Order statistics in image processing

615

2m+ 1 /~L -- m + 1

(7)

hence for m large, the mean streak length #L ~ 2. The variance a 2 of streak length is not so simple to express, however, it is easily shown that for m large, a 2 ~ 2 also (Bovik 1987). At first glance, this appears counter to the obvious fact that larger windows create extensive streaks and blotches. However, in practice, the streaking effect is not so much one where large windows create very long streaks, but instead, larger windows give rise to a larger number of moderate-length streaks which often occur in succession. For two-dimensional median filter geometries, the analysis is much more complex owing to the varying degrees of overlap between filter windows that are positioned at different orientations and distances from one another in the image. Ideally, one would like to be able to express the average "size" of a blotch, but the complexity of the problem, as well as the fact that different window geometries will yield different results has so far left this unaccomplished. In a limited analysis (Bovik 1987), the probability that two overlapping windows would yield identical outputs was calculated for a variety of window geometries. Supposing that y, and YP are the outputs of a two-dimensional median filter at coordinates n and p, where the window geometry spans N = 2m + 1 pixels, and where the windows centered at n and p overlap over exactly K _ < N pixels, then defining Sl = max(1,m + 2 - K), s2 =- min(n + 1,2m + 2 - K), the probability that the median outputs will be equal (hence both potentially part of a blotch) is

K ~(2m+sl-K)2 (m-~+~) 1(-1

Pr{yn = yp} -- 4m + 2 - K

- 1

s=sl

(4m+l-K~

(8)

\ m+s-1 J

This expression does not reduce easily, but Bovik (1987) tabulates these blotch statistics for three types of median filter window geometries (SQUARE, CROSS, XSHAPE) shown in Fig. 1, and for window spans ranging from 3 to 15 (pixel coverage depending on the window). It is found, not unexpectedly, that the CROSS and X S H A P E geometries tend towards strong horizontal/vertical streaks and strong diagonal streaks, respectively, where the S Q U A R E geometry yields streaks with little directional preference, hence tends instead to create "blotches". This blotching effect may be observed in several of the examples in Fig. 3. Refer to (Bovik 1987) for a series of image processing examples illustrating this effect. The effects of streaking and blotching are undesirable, since they may lead an observer (or an automated processing system) to believe that there is an image structure where there is not any physical correlate. Image streaks and blotches can be reduced somewhat by postprocessing the median-filtered image with a short duration low-pass linear filter (Rabiner et al., 1975), or by instead, defining the filter output to be a linear combination of the windowed signal samples, rather than just the median (Bovik et al., 1983; Bednar and Watt 1984; Lee and Kassam 1985). The latter of these defines the so-called OS filters, to be defined shortly.

616

s. T. Acton and A. C. Bovik

Output distributions of 2-D median filters As is probably evident from the previous analysis of median filter streaking, the analysis of the second-order statistical properties of two-dimensionally medianfiltered images is quite challenging. Nevertheless, Liao, Nodes and Gallagher (1985) compute expressions for the multivariate output distributions of the separable median filter, and two-dimensional SQUARE median filters. Since the expressions obtained are truly cumbersome, the reader is referred to this reference, Computation of the moving median Another area of concern regarding the median (and other OS-related filters) is the computational cost. As the window size N increases, the sorting process increases in execution time as N! and becomes prohibitively slow on a standard serial computer. The expense of arithmetic ranking has spurred the development of fast algorithms to accelerate the computation (Huang et al., 1979; Ataman et al., 1980) and special-purpose architectures tailored to the median filter operation (Oflazer 1983). The algorithm due to Huang et al. exploits the fact that many of the pixels are shared between adjacent or overlapping filter positions, and so need not be ranked at each filter position. Other fast implementations utilize the threshold decomposition property of the median filter (discussed in another chapter in this volume), which for quantized data, allows all operations to be expressed as binary median filtering (i.e., binary majority filtering).

Rank- Order filters With the broad acceptance of the median filter as an effective image smoothing device, modifications and generalizations began to appear. The simplest extensions are the rank-order (RO)filters, which deliver the k th order statistic at the filter output. Nodes and Gallagher (1982) studied the basic time-domain properties of the RO. Using the notation for the median filter, the output of the k th rank-order filter is defined as follows: if {Yn} = rank~{xn}, then y, = x(k):, .

(9)

where 1 < k < N. RO filters of interest include the max filter and the min filter. The max filter, also called the dilation filter since it dilates signal peaks, is defined by k = 2m + 1. The rain filter, also called the erosion filter since it erodes the peaks in a signal, is implemented with k = 1, the first or lowest ranking OS. The root signals produced by 1-D RO filters are described in (Eberly et al., 1991). RO filters other than the max, median, and rain find relatively practical applications. In a later section, the erosion and dilation filters will be more extensively discussed; it will be found that erosion and dilation operations can be combined in a variety of ways to create smoothing filters having interesting properties. In fact, an entire discipline which we may term digital mathematical morphology has arisen which is largely based on these two simple operations.

Order statistics in image processing

617

Detail-Preserving filters An important ingredient in assessing the efficacy of a median filter is the effect that the geometry of the window used has on information-bearing image structures. Figure 4 is instructive on this point. Figs. 4(a) and 4(b) depict close-ups of small patches of an image that are to be median filtered (for simplicity, the subimages shown are bi-valued, but the observations made here are generally applicable). The curvilinear structures that appear in the image patches could represent significant image detail, such as parts of text, minute detail in a fabric, etc. Reducing the clarity of these features as part of a noise-reducing process could greatly deteriorate interpretation of the image content. Applying a SQUARE median filter of any size will completely eradicate the structure shown in both Figs. 4(a) and 4(b), leaving nothing but blank patches the filter will regard the local patterns as "outlying values." However, applying a CROSS median filter of moderate size will leave the pattern in Fig. 4(a) completely untouched, while an XSHAPE median filter will leave the structures in Fig. 4(b) nearly unscathed. However, if the filters are reversed, viz., a CROSS median filter applied to 4(b) or an XSHAPE filter applied to 4(a), the patterns would again be destroyed. Clearly, it would be desirable to be able to incorporate selection of the window geometry into a design procedure for median filters. Apart from the difficulty of finding such an optimal geometry for a given image (according to some criterion, which would require identifying the "desirable" structures to be preserved), the problem is made more difficult by the fact that images are generally nonstationary, and apt to contain different types of features at different spatial locations. With these problems in mind, several authors have proposed "detail-preserving filters," which possess qualities of the median filter (and usually are defined in terms of combined median operations), but with additional sensitivity to structures (of variable orientation) such as depicted in Fig. 4. Filters in this "class" include multistage median filters, which combine the responses of several

m

m

IIm

mum

II mml Ili i Ilmm i

|

(a)

(b)

Fig. 4. Depictionof two simpleimage patterns for whichmedianfilteringwith differentwindowsyields radically differentresults.

S. T. Acton and A. C. Bovik

618

directionally-oriented median filters at each coordinate (Arce and Foster 1989); FIR-median hybrid filters, which combine linear filters with median filters (Neiminen, Heinonen and Nuevo 1987); and max/median filters, which combine median operations with other RO operations (Arce and McLaughlin 1987). Because of the similarity in performance of these filters (Arce and Foster (1989) perform a careful statistical analysis), only a few members of the class are described here. The class of filters can be divided into two types, the so-called unidirectional multistage filters and the bidirectional multistage filters. The unidirectional multistage filters utilize four one-dimensional median filters outputs, ,which are then combined. Fig. 5 depicts four one-dimensional windows, labeled Bi through B4. These include the ROW and C O L U M N windows from before, plus two diagonal windows. The unidirectional filters are defined as follows. At each coordinate n, denote the output of a median filter using window B i a s Yn,i for i = 1,2, 3, 4. Denote the maximum and minimum of these outputs, respectively, by

Yn,max =max{yn,i;i= 1,2,3,4} and Yn,min = min{yn,i;i = 1,2, 3, 4} .

(10)

The unidirectional multistage max/median filter is then defined to be the median of these and the current input (over which the windows are centered): Y.,max/mod = med{x.,y.,m,×,yn,mi.} .

(11)

The unidirectional multistage median filter is computed in several stages of computing medians, instead of maxima and minima. Using the same windows, define Yn,(I,3) = med{xn,yn,j,yn,3}

and

Yn,(2,4) = med{xn,Yn,2,Yn,4}

(12)

The output of the unidirectional multistage median filter is then Yn,multi-med = med{xn, Y.,(I,3),Yn,(2,4)} •

B1

B2

(I 3)

B3

B4

Fig. 5. Subwindowsused in the definitionof the unidirectionalmultistage max/median filter.

Order statistics in image processing

619

The bidirectional versions of these filters are defined simply by using CROSS and XSHAPE two-dimensional windows instead of combinations of one-dimensional windows. The bidirectional multistage max/median filter is defined as the median of three values: the current sample, and the maximum and the minimum of the CROSS and X S H A P E median filter outputs. The bidirectional multistage median filter is defined simply as the three-sample median of the current output with the two filter outputs. Arce and Foster (1989) perform an extensive statistical and empirical analysis comparing the various proposed detail-preserving filters. In addition to computing breakdown probabilities of each of the filters, empirical mean-square error and mean absolute error figures are computed from images of real scenes as well as test images. They found that the relative performances of the various filters studied varies as a function of the level of signal vs. noise (SNR), the window spans, and the noise type. Generally it was found that the detail preserving filters have considerable advantage relative to the median filter for preserving image detail, unless the noise level is very high (when the detail is obscured). Tables and detailed comparisons can be found in (Arce and Foster 1989).

OS filters The so-called order statistic filters (OS filters), or occasionally L-filters, since they are really moving L-estimators, output a linear combination of the windowed order statistics. Thus, the filter definition requires the specification of a length-N vector of filter coefficients a. In image processing applications the condition 2m+l

a~ = 1

(14)

k=l

is usually applied; from a statistical perspective, this amounts to a unbiased mean condition. In simpler terms, it means that the average level of an image will remain unchanged when filtering with an OS filter with coefficients a satisfying (15). The output of an OS filter with coefficient vector a is then defined: if {Yn} = O S a { x n ) , then N

Y" = Z

akx(k):,

(15)

k=l

From an image processing perspective, the introduction of a filter coefficient set a allows the possibility of a design procedure (whereas the median, RO, and detailpreserving filters are "automatic") based on knowledge of the image and/or noise statistics. In addition to the median and rank-order filters discussed previously, several useful fall in the class of OS filters. Of course, all of these are moving versions of standard (possibly optimal) L-estimators, and so, enjoy statistical behavior partially described by the static properties that are amply available in the statistics literature. For example, the midrange filter is implemented with a~ = 1/2 for k = 1, a~ = 1/2 for k = 2m + 1, and a~ = 0 otherwise. Standing

620

S. T. Acton and A. C. Bovik

apart as the only linear OS filter, the average filter, is constructed using ak = 1/(2m + 1). An L-inner mean filter (a moving trimmed mean, in essence) may be defined by setting a~ = 1/(2L + 1) for (m + 1) _< k < (m + 1 + L ) and ak = 0 otherwise. Fig. 6 depicts an OS filter application. The image in Fig. 6(a) was immersed in additive Laplacian-distributed noise of standard deviation a = 20.0, thus yielding Fig. 6(b). An OS filter was designed such that the weighting vector a is triangular in form. This device, referred to as the A-OS filter (read "triangle OS filter"),

(a)

(b)

(c)

(d)

Fig. 6. (a) Original "'Old Central" Image. (b) "Old Central" corrupted by Laplacian-distributed additive noise (a -- 20.0). (c) 5 x 5 SQUARE A-OS applied to noisy image. (d) 3 x 3 SQUARE W M M R - M E D filter applied to noisy image.

Order statistics in image processing

621

emphasizes OS that neighbor the window median intensity, while suppressing the extreme windowed values. The result of A-OS filtering the image using a 5 x 5 S Q U A R E window is depicted in Fig. 6(c). Clearly, the noise content of the image has been effectively removed; however, there is some expense of image blurring. Optimal O S f i h e r s

If the image samples are corrupted by additive i.i.d, noise samples, and if the image is considered sufficiently smooth that it can be modeled as nearly constant within each window position, then standard techniques may be used to design the OS filter under, for example, a standard least-squares criteria (Lloyd 1952; Sarhan 1954, 1955a,b) or using some statistical robustness criteria (David 1982; Crow and Siddiqui, 1967; Gastwirth and Cohen 1970). Such topics have been considered in (Bovik et al., 1983; Restrepo and Bovik 1988). Hence, for a sufficiently smooth signal (viz., approximately constant within the span of a window) immersed in an additive (or multiplicative) i.i.d, noise process, it is possible to construct an optimal OS filter under the least-squares or MSE criterion, and such a filter can be used with great effect to smooth signals, without the streaking effect introduced by median filtering (Bovik et al., 1983). It is worth noting that since the (linear) average filter is also an OS filter, it is always possible to derive an OS filter (under these conditions) that produces a MSE no larger than that of the optimal linear filter! When the noise becomes non-independent, or when the signal becomes fastvarying, then this simple procedure breaks down. With a little effort, however, it is possible to compute an optimal MSE OS filter to estimate a non-constant image immersed in i.i.d, noise. However, the computation involved to calculate the optimal coefficients is non-trivial. The following discussion is based on the work of Naaman and Bovik (1991). Here, an observed image signal {x.} is defined by the pointwise sum xn = S n q- r/n, where {s.} is the uncorrupted signal and {~/n} represents a constant-mean noise sequence. The mean-squared error (MSE) between the OS-filtered signal y. and the signal sn is e. = E I(yn- s.) z] = a V H n a - 2sna v/,. + s 2

(16)

Hn=E[x(n)X/n)] ,

(17)

where

is the O S correlation matrix. Its elements are given by Hi,j:n = E[x(i):nX(j):n] .

(is)

The O S mean vector is defined -- Six/.)]

,

(19)

S. T. Acton and A. C. Bovik

622

with elements /&n = E[x(i):n] .

(20)

T o create an optimal OS filter for this image processing problem, the M S E is minimized over the OS coefficients a. F o r an image containing K pixels, the M S E can be minimized over the entire set of signal samples, or over a subset o f the samples. F o r a generalized solution over K components, the M S E is K

e= Z

e.~ = a r H a

- 2a~l~ls + s r

(21)

k=l

where K

= Z ~I"k (NxN) ,

(22)

k=l

1VI=

[/~nl

I~.=l.-I~.k]

(NxK) ,

(23)

and s = (sn., &2, • • -, s.k) r

(24)

In this paradigm, the signal c o m p o n e n t s to be estimated, s, m a y or m a y not contain shared (overlapping) estimates. So, the signal m a y be estimated from K overlapping windows given by the vectors Xnk = s.~ + !1. k

(25)

x. k = (Xl:.k,Xe,k,... ,Xx:.k) r ,

(26)

S.~ = (Sl:.k, S2:nk,.'',SN:.k) r ,

(27)

!1.~ = 011:.~,t/2:.k,...,qN:nk) r

(28)

where

and

A l t h o u g h the estimation process itself is not affected by the use o f overlapping windowed components, the optimal coefficients will be different than in the nonoverlapping case. Since all of the involved correlation matrices are positive definite, the m i n i m u m M S E solution can be realized by a straightforward quadratic optimization process (with a global m i n i m u m ) as follows: O__e_e= 2~la - 21~Is = 0 Oa

(29)

Order statistics in image processing

623

which yields the optimal OS coefficients a*: a*

= H-1Ms

.

(30)

The optimal coefficients produce a MSE of: 8min

= sT ['I -- l~i2I- 1]~I] s

(31)

where I is the K x K identity matrix. For a given i.i.d, noise process with distribution F(.) known a priori, closedform expressions for the OS correlation matrix and the OS mean vector are obtainable, but at a significant computational cost (Naaman and Bovik 1991). Direct computation of these marginal and joint statistics by numerical integration involves O(K) integrals with O(N !)-term integrands. Hence, for large values of K, as encountered in image processing, and large sizes (N) of the filter window, direct computation is not practical, The composite OS moment matrices can be estimated in polynomial time (Boncelet 1987). Using a recursive, discrete approximation algorithm, the OS filter coefficients for image processing applications can be computed in a reasonable period of time. The multivariate distribution G of m order statistics x(k~),...,x(km) corresponding to a sample of independent random variables x = ( x l , . . . ,xu) r of size N _> m can be estimated using G ( k l , . . . , km; N) = Pr{x(i) < ti,

i = k l , . . . , km

m+l

= ZG(ki,...,kj-,,kj-

1,...,km-1;tl,...,tm

N-l)

j=l

× Pr{tj-l
tj}

(32)

where kl < k2 < ... < km are integers and - e c = to < tl < ...
624

S. T. Acton and A. C. Bovik

A(2, a) = e + 2~(1 - are)

(33)

where e is a unit vector of length N and 2c is a Lagrange multiplier. Differentiating (33) with respect to the coefficient vector a (and equating the result to zero), the optimal coefficients are computed as: 2c- 1 a* = H~-IlVls + ~-U~- e

(34)

which combined with (14) yields 2 ( 1 - er~r: llVls) )~c =

(35)

er/~rc i e

The results can then be used to determine a new minimum MSE. At this point, methods for computing optimal OS coefficients by a standard global MSE technique and a local structure-preserving least squares technique have been described. Using these solutions, optimal OS filters for smoothing digital imagery and other signals can be defined. The OS filter has proven to be a powerful tool in image enhancement applications. The impulse rejection and edge preservation properties of OS filters cannot be matched by linear filters, since signal discontinuities lead to overlapping frequency spectra (Restrepo and Bovik 1986). OS filters are also translationinvariant and preserve linear signal trends (Lee and Kassam 1985). W M M R filters The weighted majority with minimum range (WMMR)filters are OS filters that weight the ordered values in a subset of window samples with minimum range (Longbotham and Eberly 1993). Typically, the set of m + 1 samples with minimum range is found within a filter window of 2m + 1 samples; then the m + 1 samples are combined linearly according to a weight vector a. For example, the weights may select the median (WMMR-MED) or the average (WMMR-AVE) of the m + 1 values. The W M M R concept is a generalization of the trimmed mean filters, which "trim" away a set number of the largest and smallest order statistics in the window and then average the remaining samples (Bednar and Watt 1984). The W M M R filters have proven to be very effective edge enhancers and smoothers for signals that are piecewise constant (PICO). For an example, see Fig. 6. The image in Fig. 6(b), corrupted with Laplacian impulse noise, has been enhanced by a 9 × 9 W M M R - M E D filter (see Fig. 6(d)). Note the very sharp retention of image features despite the excellent noise suppression. An interesting root-signal analysis has also been developed for the W M M R filters. A one-dimensional signal is defined to be piecewise constant of degree d, or PICO-d, if each sample belongs to a constant segment of at least length d. Using the W M M R filter, similar concepts to linear lowpass, highpass, and bandpass filtering for PICO signals have been conceived (Longbotham and Eberly 1992). Just as Gallagher and Wise (1981) derived a root signal theory for median filters,

Order statistics in image processing

625

Longbotham and Eberly in like fashion provide a root signal characterization for the W M M R filter. In order to achieve a PICO-(m + 1) root signal by iterative application of the W M M R filter, they show that the weights a must be nonnegative and sum to unity (normalized). Also, the first and last weights (al and am+l) must be unequal (Longbotham and Eberly, 1993). Although convergence to a PICO root signal for 2-D signals is not guaranteed, the W M M R is an extremely powerful tool for smoothing images of natural and man-made scenes that consist of approximately piecewise constant regions. The W M M R filter has been successfully applied to problems in biomedicine (extracting dc drift) and in electronic circuit inspection (removing impulse noise from digital images of printed circuit boards) (Longbotham and Eberly 1992).

Relationship of OS fihers to linear FIR filters One may note that the definition of the OS filter is very similar to that of afinite impulse response (FIR) linear filter. An FIR filter is simply a linear filter that is not recursive, i.e., the current output does not depend in any way on other outputs. Using our prior terminology, a one-dimensional linear FIR filter with associated coefficient vector b has an output defined as follows: if {z.} = FIRb{x,}, then N

Zn = Z akxk:n k-1

(36)

Hence, the signal samples are not placed in rank order prior to weighting. Thus, an FIR filter weights the signal samples within a window according to either spatial or temporal order, whereas the OS filter employs an additional algebraic sorting step so that the samples are weighted according to their arithmetic rank. The hybridization of FIR and OS filters, such as the FIR-median filter (Heinonen and Nuevo 1987), is a current topic of research. The difference or similarity in performance between FIR and OS filters is signal-dependent. For a signal that is sufficiently smooth, in the sense of being locally monotonic (LOMO), OS and FIR smoothing filters operate identically. This result is formalized in Theorem 3 for 1-D filters: THEOREM 3. Longbotham and Bovik (1989): The outputs of a linear FIR filter and an OS filter are equal: OSa{xn} -- FIR,{xn}

(37)

for all even-symmetric coefficient vectors a of length N = 2m + 1, if and only if (x.} is LOMO-(2m + 1). This Theorem, which is fairly obvious, makes it clear that the description of signal smoothness in terms of local monotonicity has significance for not only the median filter, but for all OS filters. While setting the foundation for analytical

626

S. T. Acton and A. C. Bovik

characterization of OS filters, these results have also deepened the significance of the property of local monotonicity as a measure of signal smoothness. A more general and provocative Theorem is given next: THEOREM 4. Longbotham and Bovik (1989): Suppose that the elements of a length N = 2m + 1 even-symmetric coefficient vectors a satisfy a I = a2 . . . . . a2m+l-k, where m + 1 < k < 2m and one of the following: (a) al < min{a2m+2-k, • • •, am+l } (b) al > min{a2m+2_/~, • . ., am+l }. If {x,} contains at least one monotonic segment of length k, then: OSa{xn} = FIRa{xn}

(38)

if and only if {xn} is LOMO-(k + 1). An OS/FIR filter pair satisfying Theorem 4 are said to be LOMO-(k + 1) equivalent. Some interesting filters can be immediately seen to be LOMO-equiv-

alent: COROLLARY. Longbotham and Bovik (1989): Linear averaging filters of length 2k + 1 and k-inner mean filters are LOMO-(k + 1) equivalent. Many signal processing operations are simple to define in one dimension, but are difficult to extend to 2-D or higher dimensions, such as Kalman filtering. Most digital filtering operations have disparate definitions in 1-D and in 2-D (or M-D), since the filter operation must be redefined relative to the windowing geometry or other spatial criteria. This is not true with OS filters. Once each element is ranked, the OS filters are defined in the same manner for the one-dimensional, twodimensional, and M-dimensional cases; they remain as conceptually simple (although the theory may not extend across dimensions!). As a result, OS filters may be discussed in terms of general signal processing; the one-dimensional, twodimensional, and M-dimensional cases do not have to be treated separately.

3. Spatial/temporal extensions Despite the broad success of OS-retated filters in the area of image processing, which resides primarily in the fact that OS filters provide image smoothing performance that linear filters cannot, one must be careful in applying the OS technique blindly. For example, an OS filter does not capitalize on the original (spatial or temporal) ordering of the signal or image within a filter window. In this respect, there are certain processing tasks that are not well-suited for a standard OS filter. The idea behind several of the OS generalizations has been to combine the spatial/temporal ordering information of linear filters with the rank order information contained in the responses of OS filters.

Order statistics in imageprocessing

627

C, LI and permutation filters The so-called combination filters (C filters) (Ghandi and Kassam 1991) and Ll filters (Palmieri and Boncelet 1989), attempt to fuse OS filters (denoted by L) with linear filters (denoted by I). The C and Ll filters use a weighted sum of the input samples to compute the filter output. The weights are defined to be a function of both rank and temporal/spatial position. The results that have been reported indicate that improved signal estimation performance may be gained relative to the standard OS approach (i.e., the median filter) for removing signal outliers. Unfortunately, the C and L1 filters may perform poorly, or be difficult to design, for more complex or nonstationary signals, where signal discontinuities such as image region boundaries exist (Barner and Arce 1994). Consequently, the relationship between spatial/temporal order and rank-order remains a current topic of research within signal processing. In the same spirit as C and Ll filters, permutation filters (P filters) generalize the OS filter by fostering the utilization of the spatial or temporal context of a signal (Barner and Arce 1994). Since the mapping of a spatially or temporally ordered vector to a rank-ordered vector is essentially a permutation, permutation theory is used to develop the filtering theory. Each permutation filter output is defined to be the order statistic associated with mapping the original vector to the ordered vector. The P filters do not use weighted sums; each filter output is restricted to lie within the set of input samples. The operation that maps the spatially/temporally ordered vector x, to a rank-ordered vector x(,) is called the observation permutation, Px. The sample given as the output of the P filter is a function of p,,. If SN is the set of possible permutations, and H is a set whose elements (called blocks) are possibly empty subsets of SN, then the P filter is defined as follows: Fp(x.; H) = x(tx) ,

(39)

where Px E Hi, (Ix is the index of the block that contains px). So, in this construction, specific outputs are associated with permutations in the blocks of H. The arrangement of H can be tailored to specific applications such as estimation. The initial results for these very recent filters show that P filters can track discontinuities of nonstationary signals and reduce outliers effectively while also allowing effective frequency selection. The complexity of P filters is an obvious drawback, as the number of possible permutations in SN explodes combinatorially with increased window size. Reduced set permutation filters (RP filters) are P filters which consider several permutations as isomorphic and therefore equivalent, thus reducing the number of possible instances. Rank-order filters, weighted rank-order filters, and stack filters have a limited permutation transformation set, as compared to the general class of P filters; in fact they are a subset of the P class (Barner and Arce 1994).

628

S. T. Acton and A. C. Bovik

Stack filters Wendt, Coyle, and Gallagher (1986) developed the so-called stack filters. The power of these filters lies in their versatility and efficient implementation. They form a class of filters which include the RO filters as well as many of the OS filters. A VLSI implementation of the stack filter class is facilitated by exploiting the stacking property and the superposition property, also called threshold decomposition. It is shown in (Fitch et al., 1985) that RO operations can be implemented by first decomposing a K-valued discrete signal (by thresholding) into K binary signals, then applying binary filter operations (such as the binary median) on each of the binary signals, and finally "stacking" the outputs of the Boolean operations to form a K-valued output signal. Filters that can be implemented by threshold decomposition and stacking belong to the class of stack filters. The conditions under which a stack filter preserves edges and LOMO regions are given in (Wendt et al., 1986). In addition, the rate of convergence and the basic statistical properties of the stack filter are provided. Since stack filters are easily realized on a special-purpose chip and incorporated into a real-time signal processing system, research into the stack filter theory has been quite extensive. More than that, optimal stack filters can be effectively and very naturally designed according to the mean absolute error (MAE) criterion (Coyle and Lin, 1988); indeed, optimal stack filtering under the MAE criterion bears analogies to optimal linear filtering under the MSE criterion: both procedures are made possible by superposition properties, and both have tractable design procedures. Indeed, stack filtering would almost certainly find nearly ubiquitous application, if it were not for the fact that the threshold decomposition does not easily yield to intuitive or simple design. Nevertheless, impressive results have been obtained in image filtering applications (Coyle et al., 1989). The other main drawback of the stack filters are their constrained structure; these limitations are discussed in (Barner et al., 1992). Recent extensions to the stack filter theory include the creation of signal-dependent adaptive stack filters (Lin et al., 1990).

4. Morphological filters Although the connection is not always made, OS filters are also closely related to a relatively new area of image processing called morphology. As the name implies, morphological operations alter the "'shape" of objects in signals and images and are particularly useful in the representation and description of image features such as regions, boundaries, and skeletons. In a fundamental paper unifying the theory between OS filters and morphological filters, Maragos and Schafer (1987) demonstrated that any OS filter can be implemented via morphological operations. It will be convenient to introduce a new notation for the basic morphological operators. The morphological operations are highly dependent upon the size and shape of the structuring element, or window, B. Given the structuring element B, if k = 1, then {yn} = rankk{x.} = {x,}®B, the so-called erosion of {Yn} by B.

Order statistics in image processing

629

Likewise, if k = N = 2m + 1, then {y.} = rankk{x.} = {Xn} ® B which is the dilation of {x.} by B. The erosion operation can be viewed as finding the infima of the windowed sets; likewise, the dilation operation consists of computing suprema of the windowed sets defined by the structuring element. It is noteworthy that erosion distributes over arbitrary infima and dilation distributes over arbitrary suprema. In this respect, dilation (or erosion) can be implemented by the Minkowski addition (or subtraction) of the image {Xn} and the structuring element B. For binary-valued images, the way in which dilation affects the foreground objects of an image, erosion similarly affects the background. This duality relationship may be summarized by the following theorem and corollary, where X = {x,}, Xc is the logical complement of X, and X T is the transpose of X: THEOREM 5. (Haralick and Shapiro 1992): (X®B) c = Xc ® B T. COROLLARY. (Haralick and Shapiro 1992): (X ® B) c = XC®BT. For gray-scale images, the erosion and dilation operations have the dual effects of eliminating positive-going noise impulses and negative-going noise impulses, respectively. However, each also biases the signal either upward or downwards. An example illustrating this point is shown in Fig. 7. The dilated image of Fig. 7(b) has highlighted, expanded bright regions, while the eroded image (Fig. 7(c)) has been biased negatively, producing a dark, lower contrast image. Bias-reduced operators can be defined by concatenating opposite operations. The open operation is defined by the dilation of the erosion as follows: {x.} o B = ({x.}®B)

® B .

(40)

The opening is (Heijmans 1994b): (1) (2) (3) (4)

increasing such that if {Xn} C_ {y.}, then {x.} o B _C {y.} o B; translation invariant; idempotent, ({Xn} o B) o B = {Xn} o B; anti-extensive, {x.} o B c_ {x.}.

An open operation performed on an image will smooth the image, preserve edge information, and reject positive-going impulses. The closing is defined by {Xn}'B = ({xn} @ B)®B .

(41)

The close operation possesses the same four properties as the open operation, except that it is e x t e n s i v e rather than anti-extensive; therefore, {Xn} C_ {x.} • B. Note that dilation and erosion are also increasing and translation invariant. When applied to an image, the close operation will again smooth noise without removing edges but will eradicate negative-going, not positive-going, impulses. This characteristic difference between open and close can be exploited in the

630

S. T. Acton and A. C. Bovik

(a)

(b)

(c) Fig. 7. (a) Original "'Truck" Image. (b) Dilation with 3 × 3 SQUARE. (c) Erosion with 3 × 3 SQUARE. design of simple peak and valley detectors for digital imagery. Both open and close have important morphological implications in processing and interpreting images. For example, the opening will tend to remove small image regions and will separate loosely connected objects in the image. The closing, on the other hand, will remove small "holes" and gaps in the image. Because they are biasreduced operators, they do not significantly affect image region areas. So, like other OS filters, these filters, constructed from a succession of OS filters, can be used to smooth image data without altering the image structure. However, the open and close operations by themselves have limitations, as depicted in Fig. 8. The image of Fig. 8(a) was subjected to salt and pepper noise, where 5% of the

631

Order statistics in image processing

(a)

(b)

(c)

(d) Fig. 8.

pixels have been randomly changed to white or black, simulating transmission errors, as depicted in Fig. 8(b). The close operation result (Fig. 8(c)) and the open operation result (Fig. 8(d)) contain impulse noise artifacts. This situation can be mitigated by combining open and close into further-combined operators. By using further concatenation of operations we obtain the two-sided operators: ({x.} • B) o B

(42)

({x,,} o B) • B ,

(43)

and

632

S. T. Acton and A. C. Bovik

(e)

(f)

(g) Fig. 8. (a) Original "Friends" Image. (b) hnage corrupted by 5% salt-and-pepper noise. (c) Closing with 3 x 3 SQUARE. (d) Opening with 3 x 3 SQUARE. (e) Close-open with 3 x 3 SQUARE. (f) Open-close with 3 x 3 SQUARE. (g) Median filtered with 3 x 3 SQUARE. called open-close and close-open, respectively. Both of these operations have the ability to smooth noise, especially impulse noise, of both the positive-going and negative-going type, with little bias. Morphologically, open-close tends to link neighboring objects, where close-open tends to link neighboring holes. Here, it is assumed that the objects have higher mean intensities than that of the holes. Where open and close alone failed (see Figs. 8(c) 8(d)), the close-open (Fig. 8(e)) and open-close (Fig. 8(f)) filters successively eliminated the majority of the image noise. However, the performance of open-close and close-open did not match the results given by a median filter (Fig. 8(g)), using the same structuring element.

Order statistics in image processing

633

However, the open-close and close-open operations represent the result of four consecutive processing iterations, and so represent a more intense smoothing of the image. Since morphological operators can be implemented without computing the complete set of order statistics for each windowed set (only m a x and min are needed), open-close and close-open offer affordable alternatives to median smoothing for real-time processing. Furthermore, morphological operations can be realized using simple logic circuits and can be implemented on high-speed locally interconnected parallel processors. With this in mind, Haralick and Shapiro (1992) give a morphological approximation to the median. For input image {Xn}, let {y.} = ({x.} o B), and {zn} = ({x.} • B). Each sample n of the output image {qn}, is given by qn = yn if [yn _> x.I and qn = zn otherwise . If {x.} is monotone at n, the median filter response {r.} = med{x.} at sample n is Xn = Yn = Zn = qn- If {Xn} is not monotone at n and xn is a neighborhood minimum, then Yn -< rn which implies that z. _> r.. At this point qn = Zn, so that the minimum points are changed to pixel values greater than or equal to the median response. However, since x~ is a neighborhood minimum, the new value is less than the neighborhood maximum. In this way, negative-going outliers can be eliminated by this expeditious approximation of the median. The same argument is true in the case of positive-going neighborhood maxima, where the closing q. = y. is the appropriate response. In addition to filtering, these OS-based operations are extremely effective tools that can be utilized to describe the shape of particular image regions. A morphological shape decomposition (Pitas and Venetsanopoulos 1990) gives a coarseto-fine description of a given image, similar to, but more meaningful than a succession of linear filtering. The morphological representation preserves edge localization, while linear filtering tends to blur and distort image objects. A simple shape decomposition can be created by a succession of open-close filters with increasing structuring element size (Vincent 1994). This progression is often called the alternating sequential filter (ASF) method (Serra and Vincent 1992). The Pitas-Venetsanopoulos decomposition has also been applied to several imaging applications (Pitas and Venetsanopoulos 1994). Given an image object defined by a set of samples {y.} and a convex structuring element B, the largest constant rl (called the radius) is found such that {y~} o riB ¢ 0. The first-order approximation in the shape decomposition is defined as {y(1).} = {Yn} o riB .

(44)

Then, the largest radius r2 is computed such that ({Yn} {Y(1)n}) o r2B ¢ Q, which leads to the second-order approximation {y(2).} = {y(1)n } U ({y.} {y(1).}) o rzB . Given {y(0)n } = ~ , the following recursion defines the decomposition:

(45)

634

S. T. Acton and A. C. Bovik

{y(t + 1).} = {y(t).}

u

({y.}\{y(t).})

o

r,+,B

(46)

where t _> 0 and rt is the radius of the maximal structuring element rt+lB that can - 1)n}. Using (46), a union of disjoint shape compobe inscribed in { y , , ) \ { y ( t nents is created and may be used for a high-level image analysis. The multi-scale representation can then be utilized for image segmentation, edge detection, image coding, or for correspondence processes such as stereo vision, motion analysis, and feature registration. An example morphological image application is presented in Fig. 9. Here, the morphological operators are used in an optical character recognition process. First, the gray-scale image of Fig. 9(a) is thresholded (t = 150) such that the

~... (a)

(b)

",,,

/

--'--:"

i

• ", ... /..

:..

.-"

~

(d)

(c) Fig. 9.

Order statistics in imageprocessing

t'"

635

(e)

Fig. 9. (a) Original "Letters" Image. (b) Binarysegmentation.(c) Close-openwith 9 x 9 SQUARE. (d) Medial Axis Transformation.(e) Inverse Medial Axis Transform. characters are black and the background white, yielding Fig. 9(b). However, an optimal threshold does not exist to provide a complete separation of the characters and the background. To remove the unwanted holes and small regions, a 9 x 9 close-open filter is successfully applied (Fig. 9(c)). Moreover, the medial axis transform (MAT) (Blum 1973; Serra 1988) is used to provide a compact representation of this binary image (see Fig. 9(d)). One can exactly recreate the original via the inverse MAT. Image morphology, which can be viewed as a method of order statistics, is growing rapidly in its theory and applications. These economical operators offer powerful image smoothers and shape descriptors. In addition to the applications presented here, image morphology can be used to analyze granulometry (size distribution of objects), to create morphological skeletons, to bound derivative operators, to measure object-to-object distances, and to define a new method of sampling based on image shapes. Indeed, only a minute portion of the related current theory and applications can be represented here.

5. Related OS applications Edge detection

Since OS filters can be defined to have powerful edge-preserving capability, they can also be employed to locate image discontinuities. Edge detection is an important task in image processing, since it provides the subdivision of an image into delineated, structurally significant regions. Many higher-level image understanding or vision tasks depend on the success of edge detection. Unfortunately, many edge detection schemes are sensitive to noise and are expensive to implement.

636

S. T. Acton and A. C. Bovik

Bovik and Munson (1986) proposed an edge detection scheme that is both inexpensive and resilient to outliers. Their method uses median comparisons, instead of average comparisons, between local neighborhoods on each side of a prospective edge. Statistical and deterministic results show that the median-based edge detector is more effective than average-based detectors in certain circumstances. The median filter (and other OS filters) may be used to pre-process an image for traditional edge detection techniques. Median prefiltering can dramatically improve the performance of edge detectors in terms of increased noise suppression away from edges, and increased preservation of detail (Bovik et al., 1987).

Image enhancement and restoration Least squares methods The most common application of OS filters is in signal enhancement and restoration. Order statistics have also been incorporated into more sophisticated enhancement/restoration algorithms. Bovik et al. designed a image restoration algorithm that uses OS-based hypothesis testing to preserve edges and smooth between edges (Bovik et al., 1985). The order statistics within a moving window are used to detect an edge of minimum height. If an edge is present, the output value is given by the order-constrained least squares fit of the window samples. If an edge is not present, the average of the windowed samples is the output value. In this way, edge information is incorporated into the restoration process. The algorithm compared favorably to both the median and average filters, in terms of subjective perception and mean squared error. Locally monotonic regression Certain properties of OS filters have developed into a theory of signal transformations which are not OS filtering operations themselves, but are global characterizations of the result of OS filtering. Locally monotonic (LOMO) regression is a device for enhancing signals by computing the "closest" LOMO signal to the input (where closeness is defined by a given distance norm) (Restrepo and Bovik, 1993). Hence, the computation of a LOMO regression may be compared to finding the root signal produced by iterative application of a median filter; however, the signal computed by LOMO regression is optimal in sense of similarity to the input signal. Furthermore, LOMO regression yields maximum likelihood estimates of locally monotonic signals contaminated with additive i.i.d, noise (Restrepo and Bovik 1994). The high computational cost of LOMO regression is a drawback. As the number of samples in the input signal increase, the number of operations required to compute a LOMO regression increases exponentially (Restrepo and Bovik, 1993). Faster windowed LOMO methods for 1-D signals can be implemented at the expense of relaxing the requirements for optimality (Restrepo and Bovik, 1991). Unlike standard OS filtering, the extension to higher dimension signals for LOMO regression is a difficult, ill-posed problem. Approximate methods which allow small deviations from the characteristic property of local monotonicity are

637

Order statistics in image processing

described and applied to corrupted images in (Acton and Bovik 1993). The success of LOMO regression-based image enhancement has also led to the development of other forms of regression, including piecewise constant (PICO) regression. The result of a PICO regression is comparable to the result of applying the W M M R OS filter multiple times. The example in Fig. 10 depicts the ability of PICO regression to eliminate heavy-tailed noise within piecewise constant image regions. Notice the detail preservation of the PICO regression result in Fig. 1l(c). PICO regression has been utilized to enhance and restore inherently piecewise

(a)

(b)

(c) Fig. 10. (a) Original "South Texas" image. (b) "South Texas" corrupted by Laplacian-distributed additive noise (or = 11.0). (c) PICO regression of noisy image.

638

s. T. Acton and A. C. Bovik

constant signals (Acton and Bovik 1993) and also to segment images into P I C O regions (Acton and Bovik 1994a). A more difficult problem in image processing occurs when a blurring process is concurrent with noise degradation. For example, consider a digital image taken from a moving automobile, resulting in a linearly blurred scene. In addition, the image is subject to an additive thermal noise. To restore such an image, it is necessary to simultaneously sharpen (de-blur) and smooth (de-noise) the image, which are conflicting tasks. The straightforward application of an OS filter or locally monotonic regression would fail to deblur the image. Optimization-based algorithms which iteratively deconvolve the blurred image while enforcing local monotonicity to remove noise have been successfully applied (Acton and Bovik 1994b); this approach, which is very recent, seems very promising.

6. Conclusions The available theory and tools for processing digital imagery has been significantly enhanced by the addition of order statistic techniques. F r o m Tukey's discovery of the median filter to the current research on morphology, the brief history of order statistics in image processing reveals a rich, fundamental contribution. The OS filter theory provides a general approach to image enhancement and optimal nonlinear filter design. The basic concepts have been extended to capitalize on both spatial ordering and rank ordering. Also, the filters have formed the basis of a new area in image processing that explores higher level questions in image processing regarding object integrity and shape. Numerous practical applications to significant image processing problems have been reported. Nevertheless, several important theoretical questions are as of yet unanswered, and a plethora of ripe application areas have yet to be explored.

References Acton, S. T. and A. C. Bovik (1993). Nonlinear regression for image enhancement via generalized deterministic annealing. Proc. of the SPIE Symp. Visual Commun. Image Process. Boston, Nov. 7-12. Acton, S. T. and A. C. Bovik (1994). Segmentation using piecewise constant regression. Proc of the SPIE Symp. Visual Commun. Image Process. Chicago, September 25-28. Acton, S. T. and A. C. Bovik (1994). Piecewiseand local class models for image restoration. IEEEInt. Conf. Image Process. Austin, TX, Nov. 13-16. Arce, G. A. and M. P. McLoughlin (1987). Theoretical analysis of the max/median filter. IEEE Trans. Acoust., Speech, Signal Process. ASSP-35, 60-69. Arce, G. A. and R. E. Foster (1989). Detail-preservingranked-order based filters for image processing. IEEE Trans. Aeoust., Speech, Signal Process. ASSP-37, 83-98. Ataman, E., V. K. Aatre and K. M. Wong (1980). A fast method for real-time median filtering. IEEE Trans. Acoust., Speech, Signal Process. ASSP-28, 415420. Barner, K. E., G. R. Arce and J.-H. Lin (1992). On the performance of stack filters and vector detection in image restoration. Circuits Syst. Signal Process. 11.

Order statistics in image processing

639

Barner, K. E, and G. R. Arce (1994). Permutation filters: A class of nonlinear filters based on set permutations. IEEE Trans. Signal Process. 42, 782-798. Bednar, J, B. and T. L. Watt (1984). Alpha-trimmed means and their relationship to median filters. IEEE Trans. Acoust., Speech, Signal Process. ASSP-32, 145-153. Blum, H. (1973). A transformation for extracting new descriptors of shape. In: Wathen-Dunn, W., ed., Models for the Perception of Speech and Visual Forms. MIT Press, Cambridge. Boncelet, C. G. (1987). Algorithms to compute order statistic distributions. S I A M J. Sci. Stat. Comput. 8(9), pp. 868-876. Bovik, A. C., T. S. Huang and D. C. Munson, Jr. (1983). A generalization of median filtering using linear combinations of order statistics. IEEE Trans. Acoust., Speech, Signal Process. ASSP-31, 1342-1350. Bovik, A. C., T. S. Huang and D. C. Munson, Jr. (1985). Edge sensitive image restoration using orderconstrained least-squares methods. IEEE Trans. Acoust., Speech, Signal Process. ASSP-33, 1253-1263. Bovik, A. C. and D. C. Munson, Jr. (1986). Edge detection using median comparisons. Comput. Vision, Graphics. Image Process. 33, 377-389. Bovik, A. C. (1987). Streaking in Median Filtered Images. IEEE Trans. Acoust., Speech, Signal Process. ASSP-35, 493-503. Bovik, A. C. and A. Restrepo (Palacios) (1987). Spectral properties of moving L-estimates of independent data. J. Franklin Inst. 324, 125 137. Bovik, A. C., T. S. Huang and D. C. Munson, Jr. (1987). The effect of median filtering on edge estimation and detection. IEEE Trans. Pattern Anal. Mach. Intell. PAMI-9, 181-194. ZCoyle E. J. and J.-H. Lin (1988). Stack filters and the mean absolute error criterion. IEEE Trans. Acoust., Speech, Signal Process. ASSP-36, 1244-1254. Coyle, E. J., J.-H. Lin and M. Gabbouj (1989). Optimal stack filtering and the estimation and structural approaches to image processing. IEEE Trans. Aeoust., Speech, Signal Process. ASSP-37 (12). Crow E. L. and M. M. Siddiqui (1967). Robust estimation of location. J. Amer. Statist. Assoc. 62, 353 389. David, H. A. (1955). A note on moving ranges. Biometrika 42, 512-515. David H. A. (1982). Order Statistics. Wiley, New York. Eberly, D., H. G. Longbotham and J. Aragon (1991). Complete classification of roots of one-dimensional median and rank-order filters. IEEE Trans. Signal Process. 39, 197-200. Fitch, J. P., E. J. Coyle and N. C. Gallagher (1984). Median filtering by threshold decomposition. [EEE Trans. Acoust., Speech, Signal Process. ASSP-32, 1183-1188. Fitch, J. P., E. J. Coyle and N. C. Gallagher (1985). Threshold decomposition of multidimensional ranked-order operations. IEEE Trans. Circuits Syst. CAS-32, 445450. Frieden, B. R. (1976). A new restoring algorithm for the preferential enhancement of edge gradients. J. Opt. Soc, Amer. 66, 280-283. Gallagher, N. C. and G. L. Wise (1981). A theoretical analysis of the properties of median filters. IEEE Trans. Acoust., Speech, Signal Process. ASSP-29, 1136-1141. Gastwirth, J. L. and M. L. Cohen (1970). Small sample behavior of some robust linear estimates of location. J. Amer. Statist. Assoc. 65, 946-973. Ghandi, P. and S. A. Kassam (1991). Design and performance of combination filters. IEEE Trans. Signal Process. 39(7). Heijmans, H. J. A. M. (1994a). Construction of self-dual morphological operators and modifications of the median. IEEE Int. Conf. Image Process. Austin, TX, Nov. 13 16. Heijmans, H. J. A. M. (1994b). Mathematical morphology as a tool for shape description. In: YingLie O e t al., eds., Shape in Picture: Mathematical Description of Shape in Gray-level Images. Springer-Verlag, Berlin. Heinonen, P. and Y. Neuvo (1987). FIR-median hybrid filters. IEEE Trans. Acoust., Speech, Signal Process. ASSP-35, 145 153. Huang, T. S., G. J. Yang and G. Y, Tang (1979). A fast two-dimensional median filtering algorithm. IEEE Trans. Acoust., Speech, Signal Process. ASSP-27, 13 18.

640

S. T. Acton and A. C. Bovik

Kuhlman, F. and G. L. Wise (1981). On the second moment properties of median filtered sequences of independent data. IEEE Trans. Comrnun. COM-29, 1374-1379. Lee, Y. H. and S. A. Kassam (1985). Generalized median filtering and related nonlinear filtering techniques. IEEE Trans. Acoust., Speech, Signal Process. ASSP-33, 672-683. Liao, G.-Y., T. N. Nodes and N. C. Gallagher (1985). Output distributions of two-dimensional median filters. IEEE Trans. Acoust., Speech, Signal Process. ASSP-33, 1280-1295. Lin, J.-H., T. M. Selke and E. J. Coyle (1990). Adaptive stack filtering under the mean absolute error criterion. IEEE Trans. Acoust., Speech, Signal Process. ASSP-38(6), 938-953. Lloyd, E. H. (1952). Least-squares estimation of location and scale parameters using order statistics. Biornetrika 39, 88-95. Longbotham, H. G. and A. C. Bovik (1989). Theory of order statistic filters and their relationship to linear FIR filters. IEEE Trans. Aeoust., Speech, Signal Process'. ASSP-37, 275 287. Longbotham, H. G. and D. Eberly (1992). Statistical properties, fixed points, and decomposition with WMMR filters. J. Math. Imaging and Vision. 2, 99-116. Longbotham, H. G. and D. Eberly (1993). The WMMR filters: A class of robust edge enhancers. 1EEE Trans. Signal Process. 41, 1680 1684. Maragos, P. and R. W. Schafer (1987). Morphological filters-Part lI. IEEE Trans. Acoust., Speech, Signal Process. ASSP-35. McLoughlin, M. P. and G. A. Arce (1987). Deterministic properties of the recursive separable median filter. IEEE Trans. Acoust., Speech, Signal Process. ASSP-35, 98-106. Naaman, L. and A. C. Bovik (1991). Least squares order statistic filters for signal restoration. IEEE Trans. Circuits and Syst. 38, 244-257. Neiminen, A., P. Heinonen and Y. Nuevo (1987). A new class of detail-preserving filters for image processing. IEEE Trans. Pattern Anal. Machine Intell. PAMI-9. Nodes, T. A. and N. C. Gallagher (1982). Median filters: Some modifications and their properties. 1EEE 7?ans. Aeoust., Speech, Signal Process. ASSP-30, 739-746. Nodes, T. A. and N. C. Gallagher (1983). Two-dimensional root structures and convergence properties of the separable median filter. IEEE Trans. Aeoust., Speech, Signal Process. ASSP-31, 1350-1365. Oflazer, K. (1983). Design and implementation of a single-chip median filter. IEEE Trans. Acoust., Speech, Signal Process. ASSP-31. Palmieri, F. and C. G. Boncelet0 Jr. (1989). Ll-Filters-A new class of order statistic filters. IEEE Trans. Acoust., Speech, Signal Process. ASSP-37. Pitas, I. and A. N. Venetsanopoalos (1990). Morphological shape decomposition. IEEE Trans. Pattern Anal. and Mach. Intell. PAMI-12, 38~45. Rabiner~ L. R., M. R. Sambur and C. E. Schmidt (1975). Applications of a nonlinear smoothing algorithm to speech processing. IEEE Trans. Acoust., Speech, Signal Process. ASSP-23, 552-557. Restrepo (Palacios), A. and A. C. Bovik (1986). Spectral analysis of order statistic filters. IEEE Int7 Conf Acoust., Speech, Signal Process. Tokyo. Restrepo (Palacios), A. and A. C. Bovik (1988). Adaptive trimmed mean filters for image restoration. IEEE Trans. on Acoustics, Speech, and Signal Processing. 36, 1326-1337. Restrepo (Palacios), A. and A. C. Bovik (1991). Windowed locally monotonic regression. IEEE Int7 Conf. Acoust., Speech, Signal Process. Toronto. Restrepo (Palacios), A. and A. C. Bovik (1993). Locally monotonic regression. IEEE Trans. Signal Process. 41, 2796-2810. Restrepo (Palacios), A. and A. C. Bovik (1994). On the statistical optimality of locally monotonic regression. 1EEE Trans. Signal Process. 42. Sarhan, A. E. (1955). Estimation of the mean and standard deviation by order statistics. Ann. Math. Stat. 25, 317-328. Sarhan, A. E. (1955a). Estimation of the mean and standard deviation by order statistics Part IL Ann. Math. Stat. 26, 505-511. Sarhan, A. E. (1955b). Estimation & t h e mean and standard deviation by order statistics Part III. Ann. Math. Stat. 26, 576-592.

Order statistics in image processing

641

Serra, J. (1988) hnage Analysis and Morphology." Volume 2." The Theoretical Advances. Academic Press, London. Serra, J. and L. Vincent. (1992). An overview of morphological filtering. Circuits, Systems, and Signal Processing. 11(1), 47-108. Tukey, J. W. (1971). Exploratory Data Analysis. Addison-Wesley, Reading, MA. Tukey, J. W. (1984). Nonlinear (nonsuperimposable) methods for smoothing data. In: Conf. Rec., EASCON, 673, 1974 (also available in The Collected Works of John W. Tukey, II, Time Series: 1965-1984, D. R. Brillinger, Ed. Wadsworth, Monterey, CA.). Tyan, S. G. (1981). Median filtering: Deterministic properties. In: Two-dimensional Signal Processing: Transforms and Median Filters. T. S. Huang, ed. New York: Springer-Verlag. Velleman, P. F. (1980). Definition and comparison of robust nonlinear data smoothing algorithms. J. Am. Star. Ass., vol. 75, pp. 609-615. Vincent, L. (1994). Morphological area openings and closing for grey-scale images. In: Ying-Lie O, ed., Shape in Picture. Mathematical Description of Shape in Grey-level Images. Springer-Verlag, Berlin, 197-208. Wendt, P. D., E. J. Coyle and N. C. Gallagher, Jr. (1986). Stack filters. IEEE Trans. Acoust., Speech, Signal Process. ASSP-34, 898-911.