Adaptive rank order based filters

Signal Processing 27 (1992) 1-25 Elsevier

1

Adaptive rank order based filters Philippe Salembier Signal Processing Laboratory . Swiss Federal institute of Technology, CH-/015 Lausanne, Switzerland Received 17 January 1991 Revised 20 June 1991

Abstract. This paper is devoted to the adaptive optimization of a class of nonlinear filters referred to as rank order based filter . This class includes, as special cases, classical rank order filters and morphological filters with flat structuring element . The definition highlights their common origin based on sorting operations . It is shown how to adapt the filter mask and/or the rank parameter in order to minimize a criterion such as the mean square error (MSE) or the mean absolute error (MAE) . The proposed algorithms are similar to the LMS algorithm and one of their major advantages is their simplicitly . The algorithms derivation is achieved by replacing the input/output relation of the filter involving sorting operations, by an implicit formulation using only sign functions . Several practical examples for one-dimensional signals are described illustrating the convergence properties and the difference between MAE and MSE optimization . Finally, the algorithms are applied to various real images, illustrating the ability of adaptive nonlinear filtering in coping with textures, in removing different types of noise and in dealing with non-stationary signals .

Zusammenfassuag. Diese Arbeit ist der adaptiven Optimierung einer Klasse nichtlinearer Filter gewidmet, die als rangordnungs-basierte Filter (rank order based filter) bezeichnet werden . Diese Klasse beinhaltet als Sonderfalle die klassischen Rangordnungsfilter and morphologischen Filter mit flachen Strukturelementen . Die Definition beleuchtet ihre gemeinsame Wurzel, die auf Sortieroperationen beruht . Es wird gezeigt, wie die Filtermaske and/oder die Rangparameter adaptiert werden konnen, um ein Kriterium wie das des mittleren quadratischen Fehlers (mean square error, MSE) oder des mittleren absoluten Fehlers (mean absolute error, MAE) zu minimieren . Die vorgeschlagenen Algorithmen sind dem LMS Algorithmus ahnlich and einer ihrer Hauptvorteile ist ihre Einfachheit . Die Herleitung der Algorithmen erfolgt, indem die Eingangs- ;Ausgangsbeziehung des Filters, die Sortieroperationen umfalit, durch eine implizite Formulierung ersetzt wird, die nor Vorzeichenoperationen benotigt . Mehrere praktische Beispiele fur eindimensionale Signale werden beschrieben . urn die Konvergenzeigenschaften and die Unterschiede zwischen der MAE and MSE Optimierung zu veranschaulichen . Schlielilich werden die Algorithmen auf verschiedene echte Bilder angewendet, um die Tauglichkeit der adaptiven nichtlinearen Filterung bei der Verarbeitung von Texturen . bei der Entfernung verschiedener Typen von Rauschen and bei der Behandlung nichtstalionarer Signale aufzuzeigen .

Resume . Cet article est consacre a ('optimisation adaptativc d'une class¢ de filtres nonlineaires appeles'rank order based filter' . Cette classe inclut entre autre les filtres de rang classiques et les filtres morphologiques a element structurant plat . La definition met en evidence ('origin commune de ces filtres basee sur la notion de tri . On montre comment it est possible d'adapter le masque du filtre et/ou le parametre de rang pour minimiser on critere tel que I'erreur quadratique moyenne (MSE) on I'erreur absolue moyenne (MAE) . Les algorithmes proposes sent similaires a l'algorithme LMS et sent tres simples . Its sent obtenus en remplagant Ia relation entree/sortie du filtre impliquant une operation de tri, par une formulation implicite n'utilisant que la fonction sign . Plusieurs examples sur des signaux monodimensionnels sont decrits pour montrer les proprietes de convergence et la difference entre ('optimisation suivant le MSE et Ie MAE . Finalement, les algorithmes sent utilises avec des images reelles pour illuster le comportement du filtrage adaptatif nonlineaire en presence de textures, de different type de bruit et do nonstationnarites .

Keywords. Adaptive filtering, rank order filters, morphological filters .

0165-1684/92/$05 .00 © 1992 Elsevier Science Publishers B .V . All rights reserved

2

P. Salembier / Adaptive rank order based filters

1 . Introduction Median filters, or more generally rank order filters [8], are a class of nonlinear filters which have become very popular in signal and image processing because of their attractive properties. For example, they can suppress impulse noise or more generally small signal components while preserving very important features such as edges . This characteristic differentiates them strongly from linear filters which smooth impulse noise and blur edges . Another class of increasingly popular filters are morphological filters [6, 7, 13] . They are derived from a shape-oriented approach to signal analysis . They analyze the geometrical structure of a signal by locally comparing it with a predefined elementary form called a structuring element . Both types of filters share the same origin involving sorting of the input values . They have been widely used in image processing, pattern recognition, robot vision, etc . Most of the time, the choice of a specific filter is performed with the so called structural approach . It consists in analyzing and selecting a relevant set of signal structures, and in choosing a filter which is able to properly interact with these structures . For instance, in the field of median or stack filtering, this approach has created a great number of studies on root signals which have the property of being preserved by the filter [3, 8] . In mathematical morphology [13], the structural approach leads very often to rely on the geometrical structure of the original signal in order to choose a specific structuring element. Even if this approach has been successfully applied in many practical cases, a lot of research efforts are devoted to alternative solutions . Indeed, the structural approach often requires a fair amount of practical experience and intuition . Recently, alternative approaches have been proposed . They formulate the design as an estimation problem . The filter is optimized in order to minimize a criterion such as the mean square error (MSE) or the mean absolute error (MAE) . In a similar way to that of linear filter optimization, this approach yields two different types of solutions . The first one intends to find an optimal fixed filter . Signal Prnressing

It assumes stationarity of the various signals and it also requires a certain knowledge about the noise characteristics [2, 4, 16] . Another drawback of this approach is its computational complexity . To avoid these difficulties of fixed optimum filters, several adaptive solutions to the optimization problem have been studied . They reach the optimum point step by step and, at least partially, overcome the drawbacks mentioned above . One of the original works in this field deals with the large class of stack filters [5] . Taking into account the internal structure of the filter, it is shown how the boolean functions can be adapted in order to minimize the MAE between the filter output and a desired signal . However, the number of variables to optimize grows exponentially with the filter window size and this results in a rapid increase in computation . In [ 16], the internal structure of the stack filter is also considered but the optimization is not performed on the boolean functions themselves but on the socalled separating functions . These functions are, in fact, a polynomial representation equivalent to the positive boolean functions . By performing a linearization of the nonlinearity (step function) involved in the representation with separating functions, it is shown that any linear adaptive algorithm can be used to optimize the filter . Other studies have also been dedicated to the optimization of specific filtering structures such as hybrid median filters and LI filters [9], weighted median filters [ 12], etc . But, to the author's knowledge no study has been specifically devoted to classical rank order operators in the context of MAE or MSE optimizations, and no solutions are available for adaptive morphological filtering. This study is devoted to the adaptive optimization of a class of nonlinear filters referred to as rank order based filters . This class is built by composition of two dual rank order filters . It includes, as special cases, median and morphological (with flat structuring element) filters . Once a particular composition of the dual filters has been selected, the filter characteristics depend on a rank parameter and on a mask defining the location of the input values . We show in this paper how to adapt the rank parameter and



P . Salembier / Adaptive rank order bared filters

the filter mask in order to minimize the MAE or the MSE criterion . The proposed solutions can be applied for example to rank order filters and to morphological filters with flat structuring element . One of the interesting features of the approach described here is that it leads to very simple algorithms and also to fast convergence . The organization of the paper is as follows . In the second section, the class of rank order based filters is defined and some basic properties are given . The class definition highlights a common origin based on a sorting operation which is difficult to deal with in an optimization problem . To by-pass this difficulty, an alternative formulation of the rank order operator is proposed . Its main advantage is that it only involves sign functions which are more suitable for the adaptation . This expression is obtained using an implicit formulation of the filtering process . The third section deals with the statement of the optimization problem . Both MAE and MSE criteria are considered . In practice and with an approach similar to the least mean square algorithm [ 15], they lead to very similar formulas . The problem associated with the binary nature of the mask optimization is also discussed, and a solution is proposed . The fourth section describes in detail how the rank and the mask of a single rank order filter can be optimized . Practical examples, with one-dimensional signals, illustrating the convergence of the algorithms and the difference between MAE and MSE optimization, are presented . These results are generalized in the fifth section to more complex filters like morphological opening, closing or alternating sequential filters . Finally, the last section is devoted to the application of these techniques to image processing . Several examples are presented, showing the ability of such adaptive nonlinear filtering techniques in coping with different kind of texture, in removing noise and in dealing with non-stationary signals .

2 . Rank order based filters 2 .1 . Definition and properties

The filter class under investigation here is referred to as rank order based filter (ROBF) [10] .

3

The main motivations of using this class are twofold . First, it allows dealing within the same framework with well known nonlinear filters like median [8], rank order or morphological filters [13] . The class definition highlights their common origin based on rank ordering . Therefore, it allows the investigation of the optimization problem with the same approach . The second interest of this class is that it contains less popular filters which may be useful for specific applications . An example involving noise with asymmetrical probability density function will illustrate this point . ROBFs are based on composition of two dual classical rank order filters (ROFs) . Let us recall their definition . Let x ; be a discrete signal defined on a multidimensional space Z d. Denote by M a multidimensional mask containing N points (I MI=N, if I I is the set cardinality) and by j an index belonging to the mask M . Finally, let r be the normalized rank of the filter, 0 <, r <, 1 . Then, the output of a ROF at location i can be computed by taking all the points x, belonging to the mask centered on i, by sorting their values in ascending order and by taking the nth value (n is the nearest integer value of (N-1 )r+ 1) of the ordered list . More formally, the output of the rank order operator E,,,[ ] is given by Yt=Er .M[xl]= Rank {x1 - IJcM},

(2.1)

(N- I)r+ I

where Rank„{X} selects the nth ordered value of a set X. Note that a normalized rank r (r=(n-l)/ (N-1)) is used instead of the classical rank n which defines the position of the output value in the ordered list . This is for commodity reasons . Since the optimization studied in this paper involves a filter mask adaptation, the number N of points belonging to the mask is allowed to vary during the process . Therefore, it is more convenient to use a normalized rank parameter r going from 0 to 1 instead of the real rank going from 1 to N. ROBFs are based on combinations of this classical ROF and its dual which is obtained by taking vd . 27 . No . I . April 1992



4

P. Salem bier i Adaptive rank order based filters

the rank r' equal to I -r and the reflected mask M' of M with respect to its origin. Let us denote by D,,,,,[ ] this dual operator . Its definition is given by the following formula : Dr,M[x+] =

EI -r.M•[xi]

Rank

{x ;+ ;IjEM} .

(2 .2)

(N- 1)(1-0+1

The type of combinations which are done with $„u[ ] and D,._.1[ ] are similar to that done in mathematical morphology with the erosion and the dilation . Let us mention three examples :

where E„N [ ],< E;, M[ ] for all x i .

stands

for

E,,,,,[ ] ( Trivial combination),

Two ROFs

D,,_„[Er,M[ ]],

Four ROFs

E,,,,,[D,,M[D,,,,,[E,,M[ ]]]] .

2. In the case of one-dimensional signal and connected mask M, PROPERTY

2 .2 . Alternative definition of

ROBFs

As proposed in the introduction, our objective in this paper is to investigate the possibility of making these filters adaptive in order to optimize a criterion such as the Mean Square Error (MSE) or the Mean Absolute Error (MAE) . The adaptation approach chosen in the following relies on the use of the derivative of the filter output with respect to the parameter to optimize . It would have been convenient to have an analytical and easily derivable expression of the output as a function of the input signal . Unfortunately, since ROBFs share the same origin based on rank ordering, this is not possible . Indeed, the derivative of the sorting operation involved in (2 .1) is not easy to obtain or to manipulate . To solve this problem, we propose now an alternative formulation of the rank order operator . Denote by {xi _ ;}j ,, the set of input values inside the mask at location i, and by y, the nth ordered value (as discussed in Section 2 .1, n is equal to (N- l )r+ 1) . Let us call S the sum of signs of (x,- j - y,) for all .i. Its expression is given by

], D,,M[ ] CDr.M[ ]

VM Vr, r' such that r < r', Signal Processing

(2 .4)

Note that it is easy to find counter examples in two dimensions even with a connected mask [11] .

PROPERTY1

M[

1]

Vr, r' such that r 5 r' .

Note that depending on the value of r, some well known filters can be viewed as special cases of ROBF : - If r=0 .5, E,,,,,[ ] is the median filter . - If r=0, E,. .,,[ ] ( respectively D,,,,,[ ]) is the classical morphological erosion (respectively dilation) with a flat structuring element . In the same way, D,, M[E,, M [ ]] ( respectively Er, M[Dr,,w[ ]]) is an opening (respectively closing) . - If r=1, similar morphological filters are obtained by inverting the roles of E,,,,[ ] and D, .M[ ] . Some basic properties of ROBFs are given in [ 11 ], where it is pointed out that ROBFs belong to the more general class of stack filters [14] . This induces a certain number of properties such as : every ROBF is translation invariant, upper semicontinuous, increasing (the increasing property is also called the stacking property) and commutes with thresholding . Two ordering relations between filters are also stated . The first one is obviously derived from the definition .

Er.n,[ ]
M [x;]
The second property is less obvious and is valid in the restricted case of one-dimensional signals and connected masks . The reader is referred to [11] for the proof.

Dr,M[Er,M[ ]]
Single ROF

E,,

(2 .3)

S = Y_ sgn(xi_, - yi), j€M

(2 .5)



P. Salemhier

Adaptive rank order based filters

where the sign function is defined by

sgn(a)=

1 0 1

F(x,_j , y„ r) . As will be seen in the sequel, the

advantage of this implicit formulation is that it involves sign functions instead of sorting operations . This leads to a tractable formulation of the optimization problem.

if a>0, if a=0, if a<0 .

Suppose that the samples .x i _ 1 are all different (i .e. noisy signal) . For clarity purpose, let us first compute the sum S for three particular values of the rank parameter r : If r=0 . The output y, is simply the minimum of the set {x ; j } . Then, all terms (x i _,-y ;) are strictly positive except one which is equal to zero (remind that x;_ ; are different) . The sum S is then equal to N- 1 .

If r= 1 . The output y i is the maximum of the set {x,-j} . All terms (x ;_j - yi) are strictly negative except one which is equal to zero, and the sum S is equal to -(N- 1) . If r=0.5. The output y ; is the median value of the set {x,_,} . There are as many terms (x ;_ j -y ;) which are positive than negative . The sum S is then equal to 0 . If the dependence of the sum S with respect to the parameter r is approximated by a linear function (it is in fact a stair case function since S is quantized), the general expression of S is given by the following formula : S=

5

sgn(xi j-v1)=-(2r-1 )(N-1) . (2 .6) jE M

The difference between these two expressions of S is zero, and gives the following formulation of the rank order operator . DEFINITION. The output y ; of the rank order

operator E,,,ti,[ ] at location i is the value belonging to the set {z i j} jEM such that

3. The optimization problem

3 .1 . The adaptation criterion

As discussed in the introduction, there are two different ways of dealing with a filtering optimization problem . The first one involves a fixed filter which is optimum with respect to a certain criterion and for stationary signals . The second approach relies on adaptive techniques . It is often more simple and it allows to deal with non-stationary signals . Moreover, no a priori knowledge on the noise statistic is required . This study addresses the adaptive approach . It can be stated as follows . Given an original signal x,, a desired signal d ; and a general ROBE f_,,,, what are the values of the rank parameters r and the filter mask M which minimize a cost function C? To reach the optimum point an iterative algorithm can be used . Among the various approaches which can be considered, we have chosen a specific one, similar to the Least Mean Square (LMS) algorithm [15], because it is known to be simple and robust in the linear case . Of course the optimality sense depends on a specific criterion . In this study, both the MSE and MAE are considered . Suppose that y; is the filter output which can be considered as an estimate of a desired process d, . The optimal MSE filter minimizes the squared error :

F( .z;_,, yr, r) = Y sgn(x,-j-y,) jEM

+(2r-1)(N-1)=0 .

(3 .1)

CMSE=E[(d,-y,)2], (2.7)

In this definition, the explicit formulation of .1) has been lost, and the filter output is (2 expressed through the implicit function

where E is the mathematical expectation . The optimal MAE filter reduces the average absolute difference . This can be written as CMAE=E[sgn(d,-y,)(d;-yI)] . V .I .

27 .

N. 1 .

(3 .2) April 1992



6

P. Salembier / Adaptive rank order based filters

LMS-like algorithms reach the optimum solution by computing an estimate of the gradient with respect to the parameter to be optimized, and by correcting the current solution according to this estimate . Because of the poor quality of the estimation, the correcting term is only a very small fraction of the gradient . This acts in fact as a lowpass filtering and the convergence is obtained in a statistical sense . Let us see the difference between MSE and MAE in the context of the LMS algorithm . The gradient estimate of the criterion does not take into account the mathematical expectation involved in (3 .1) and (3 .2) . In the MSE case, this leads to the following estimate : (3 .3)

VCMSE__-2(d,-y,)Vy, . In the case of MAE, the same approach gives VCMAF" +(d,-Y,)V sgn(d, - Y,) - sgn(dr- y,)Vy, .

(3 .4)

The gradient of a sign function is proportional to a delta function 8(d, y ,) . In practice, if the signals are noisy, the filter output and the desired process are always different. This means that the first term of (3 .4) can be neglected, giving the following approximation of the gradient : VCMAE~ - sgn(d, - yr)Vy, .

(3 .5)

Comparison between (3 .3) and (3 .5) shows that, in the context of the LMS optimization and with the approximations done here, the main difference between MAE and MSE optimization lies in the sign function of the error (d,-y,) . This difference will be illustrated in the following by practical examples . Note that sometimes the expression of (3 .5) is considered as being a simplification of the LMS algorithm in the MSE case . The last difficulty in the optimization problem is linked to the discrete nature of the parameter to optimize . Concerning the rank parameter, the problem can be solved by optimizing directly the continuous normalized rank r, and defining n as being the nearest integer value of (N- 1)r+ 1 . But Signal Processing

the mask optimization is less straightforward . Indeed, the mask shape can be defined by a set of boolean or binary values which indicates whether or not the corresponding location has to be taken into account in the computation of the output .

3 .2 . From binary optimization to continuous optimization Iterative algorithms analogous to the LMS are suitable for continuous optimization problems . Indeed, they converge to the optimal solution by performing a great number of small corrections on the current estimate . Obviously, they are not directly applicable for binary parameters such as those used to define the shape of the filter mask . This problem can be solved by the following procedure : a search area, which should contain the optimum mask, is defined . A continuous variable is assigned to each possible location in the search area. Finally, the current filter mask is obtained by thresholding the set of continuous parameters by an arbitrary value which is assumed to be equal to zero in the following . If the continuous value corresponding to a particular location of the search area is greater (respectively lower) than zero, this location is considered as belonging (respectively not belonging) to the filter mask . This procedure is illustrated in Fig . I for a one-dimensional filter .

--aoao-¢00000 0

_Points belonging to the rdw mask Points not belonging m the filter mask

oDG-~Da C

Continuous values to adapt i 0

5

10 Index

'5

20

Fig . 1 . From binary optimization to continuous optimization .



7

P. Salembier / Adaptive rank order based filters

The goal of the adaptation process is now to optimize the filter mask by modifying the set of continuous values of the search area instead of the binary values of the mask . In this way, the problem becomes a continuous optimization problem . Note that even if the parameters to adapt are the set of continuous values, the optimization criterion still concerns the mask itself . Let us give now a formulation of the ROF, equivalent to (2 .7), but which takes into account the search area and the set of continuous values . Suppose that A denotes the search area and ;),} ;€A is the set of corresponding continuous m values . As in the previous sections, x, and y, represent the filter input and output, respectively . Let us compute the sum S of the signs of the differences (x,_ ;-y,) for all x,_ ; belonging to the filter mask . It is equal to the weighted sum of the signs of (x,_ ;-y,) for all x, j belonging to the search area, the weights being equal to one or zero depending on whether or not the location belongs to the filter mask . The sum S can he expressed by S= L sgn(x,_,-y,) ,oAt

_ Y

;(sgn(ml )+ 1) sgn(x,_ -y,) .

(3 .6)

When the normalized rank parameter varies from zero to one, the sum S varies from

Y '(sgn(ng)+l)-l je .4

to

-I '(sgn(ml)+l)+l . leA

As in Section 2 .2, the sum S is given, in the general case, by

S=-(2r-l)[ ,

A

i(sgn(mj )+1)-1 . J

(3 .7)

Combining (3 .6) and (3 .7) leads to the following implicit formulation of the rank order operation :

DEFINITION . The output y , of the rank order operator E, .,,,[ ] at location i is the value belonging to the set {x,_ ;},_ ., such that F(nt, x,-,,

v, . r)

-Y'(sgn(n)+l) ,,A

x [sgn(x,- ; -y,) +2r - 1 ] + I -2r = 0.

(3 .8)

Let us see now how this formula can be used to perform an optimization of a single ROF .

4. Adaptation of a single ROF : E,,,,[ ] 4.1 . Adaptation process In this section, we present the adaptation algorithm for both the mask shape and the rank parameter . Several practical examples in the case of one-dimensional signals are presented . The application of these techniques to images will be shown in another section . Let us recall the formulation of the problem . x,, y, and d, denote the filter input, the filter output and the desired signal, respectively . We want to find the best rank value and the best mask in order to minimize either the MSE or the MAE between y, and d, . As discussed in Section 3 .2, the best mask can be found by optimizing a set of real values {m;} JeA corresponding to the search area A . If and r denote the current solution, the next estimates {m t, EA and r' can be computed in the case of MSE optimization with the following formulas : 1 mk = 'Ilk +2µ(d,-y,) ~ y' cm,k

r'=r+21 (d,- Y

.) cry, . Or

VkeA,

(4 .1)

where ~ and p are two parameters controlling the convergence of the algorithm. Vel .'-7. No . 1 . April 1432



8

P. Salen,bier / Adaptive rank order based filters

In the case of MAE optimization, similar formulas can be derived :

mk=mk +2µ sgn(d-yi)

aY,

can be extracted from (3 .8) : OF

1 a sgn(m k )

amk

2

[sgn(x i k-yi)+2r - 1 ] .

am,

VkeA,

(4 .5)

amk r'=r+2 .Asgn(di -y i )

a~.

(4.2)

The only problem is to find an expression of the derivative of the filter output with respect to m k and r . As stated in Section 3 .2, the implicit formulation of the input/output relation of a single ROF can be used for this purpose . Let us start with the filter mask optimization .

4.2 . Adaptation

of

the filter mask

In order to compute the derivative of the output y, with respect to a particular parameter Mk, one can use the classical technique of implicit function derivation : the total derivative of a function F of several variables with respect to Ink is equal to the sum of the partial derivative of F with respect to all its variables multiplied by the partial derivatives of the variables with respect to Mk . Applying this technique to the function Fof (3 .8), noting that all mj for j#k, all xi _j and r (assumed to be fixed in this section) do not depend on Mk, and taking into account the fact that F() is zero, leads to

Thus, OF

=8 im, [ sgn(ri_k-y,)+2r-1] .

(4 .6)

am, For the second one, (3 .8) gives OF= ay,

-E', (sgn(m,)+1)28 is

_ , .

( 4 .7)

jEA

This last expression can be simplified. In the case where all (x i j - yi ) are different from zero except one of them (see Section 2.2) . Let us denote jo the unique index j such that x i _j =yi . Then, (4 .7) reduces to OF --2 sgn(mjo)+ 1 Oyi

4 8)

2

Since jo corresponds to the location which has given the output yi , it belongs to the current mask and sgn(mjo )=1 . Equation (4 .8) is then simply OF =-2

.

(4 .9)

ayi

Finally, combining (4 .4), (4 .6) and (4 .9) gives the derivative of yi with respect to m k :

dF(mj , x,- j , yi , r)

ay

dmk

amk OF + OF ayamk

eyi

j-0 . amk

(4 .3)

The partial derivative of y i with respect to Ink is then given by ay i OF /iaF amkk -- amk ay i

(4 .4)

The computation of the derivative of y, with respect to m k is reduced to the computation of two derivatives of F with respect to m k and Iv, . The first one Signal Prowsnng

= ~8 ("")s [ gn(x,_ - k - Y +) +2r-1 l ~

(4 .10)

The presence of 8 1 ,„, in this last formula indicates that a correction should be applied to mk only if it is equal to zero! Of course, this has no meaning for practical applications . So, one can use an approximation of this last 8 function . As the delta function comes from the derivative of the sign function, several approximating functions such as derivative of piecewise linear, arctan or sinh functions have been tried . In practice they lead to similar results . For simplicity reasons the linear piecewise approximation illustrated in Fig . 2 has been chosen .



P . Salem bier / Adaptive rank order

Original sgn( .)

9

hared filters

It is interesting to remark the simplicity of this formula which is only slightly more complex than the classical LMS algorithm for the linear filtering case . Note also that the same formula can be used for all rank order filters, in particular erosion (r=0), median filter (r=0 .5) or dilation (r=1) . Let us illustrate this adaptive algorithm with some practical examples .

sgn( .) Linearpiecewise approximation ->

First example : Convergence

Derivative: delta function

derivative

Fig. 2 . Approximation of the remaining delta function .

With this approximation which only deals with the remaining S function and not the original expression of F( ), the correction is applied to Mk values lying between two bounds . In practice, the bounds values do not have a strong influence on the algorithm results . Indeed, they prevent the m k variables to be too far from the threshold . But this remoteness is highly dependent on the convergence parameter which defines the speed evolution of the variables . For practical applications reported here, the bounds are arbitrarily fixed to -1 and 1 . The final estimate of the derivative is then c'yt =, z(sgn(xr-k - Y .)+2r - 1)

am,

for - I

< mk < I .

(4.11)

of the

algorithm

The first example shows the convergence of the algorithm in the case where the filter is used as a predictor . The input signal is a sine wave with a period equal to five time samples corrupted by a zero mean noise . The signal-to-noise ratio is around 15 dB . For each location, the filter estimates the value of the current input signal knowing only a limited number of past values, The desired signal d-is then the input signal x, itself . The search area A is defined as a window of length NA which is equal to 15 in this particular test . This adaptive prediction scheme is presented in Fig . 3 . The adaptation is done accordingly to the MSE criterion corresponding to (4 .1) . The first investigated filter is the median filter (r=0 .5) . Its initial mask is equal to the search area itself (m,=1 . 1 <,j,<15) . Figure 4 presents the time evolution of the absolute prediction error for three values of the convergence parameter p of (4 .1) . In all cases, the filter mask converges towards its optimal value which consists of the set of points

Adaptive algorithm mZ

m [Al

ROF

Fig. 3 . Adaptive prediction with a ROF . Vol . 27. No. I, April 1992



10

P. Salembier / Adaptive rank order based filters

Prediction error

Prediction error

. 1 11=0.0001

0

0

In- -

µ=0.001

µ=0 .01

0 0

2000

0

2000

Iterations

Iterations

Fig . 4 . Mask optimization in the case of median filter .

corresponding to a multiple of the sine wave period . The optimal prediction is achieved with the following rule : z

f=Median[x ;_ s, x i _ io ,

x,_, 5] .

Prediction emir

(4 .12)

An interesting point, which is clearly visible in the case of slow convergence (p =0 .0001), is that the absolute prediction error reaches its optimal value by steps . This behavior is due to the binary nature of the mask optimization problem . Each jump of the absolute prediction error corresponds to the introduction or the omission of specific locations in the mask. This contrasts with linear adaptive algorithms which produce a smooth exponential convergence of the prediction error [15] . Comparison of the three convergence rates shows that the residual error when the filter has converged, does not heavily depend on the p parameter . This also contrasts with the linear case where a high convergence speed is achieved at the expense of a rapid increase of the residual prediction error. Figures 5 and 6 present the same experiment with an erosion (r=0) and a dilation (r= 1), respectively . As can be seen, the convergence is also achieved by steps . As far as the prediction error is concerned, the performance of the algorithms is close to the one obtained with a median filter . Signal Processing

Fig. 5 . Mask optimization in the case of erosion .

µ=

0

l

0

0 .0001

µ = 0.001

0

2000 Iterations Fig. 6 . Mask optimization in the case of dilation .

Finally, Fig. 7 compares the time evolution of the absolute prediction error and of the continuous values associated to each location of the search area . In fact only the five first continuous values are displayed . Indeed, the input signal being periodic, the remaining values have the same kind of evolution as the five first ones . The filter is a ROF with a normalized rank equal to 0 .25 . In this figure, one can see that the evolution of the continuous values is rather smooth (as discussed above, they



P . Salembier / Adaptive rank order based filters

II

Signal

Coefficients m5

.... ................. ... ... .

A I'

MI

1 ,1

1112

M3

\

V J

I

MAE adaptation

J

Prediction error

J

1

1

Original

too Time 0

Fig . 8 . Comparison MSE/ MAE with a median filter . Iterations

Fig. 7 . Comparison of coefficients and error convergences .

well known that in such a case the MAE criterion is more robust than the MSE criterion . 4.3 . Adaptation

are bounded by -I and 1) . As suggested above, each zero crossing corresponds to a jump of the prediction error .

Second example : Comparison adaptation

of MAE

and MSE

The next experiment illustrates the difference between MSE adaptation (eq . (4 .1)) and MAE adaptation (eq . (4 .2)) . The test signal is the same noisy sine wave used before, but one period has been totally removed and replaced by a random signal . This original signal is presented in the lowest part of Fig . 8 . An adaptive median filter is used to restore the signal . Both MSE and MAE adaptive algorithms are applied, and the convergence parameters are set to have the same convergence speed for both algorithms . As can be seen, the restoration obtained with the MSE algorithm is slightly affected by the presence of the missing period, whereas in the case of MAE optimization the sine wave is perfectly reconstructed . This result is in accordance with intuition since the missing period can be viewed as an impulse noise, and it is

of the

rank parameter

In order to optimize the normalized rank parameter r, one can use the same approach as for the mask shape : computing the total derivative of the implicit function F with respect to r gives

ay,

OF j8F

Or

Or ;

av,

(4 .13)

The only unknown partial derivative is aF/or, which can easily be derived from (3 .8) : OF =

Or

2 Y4i(sgn(m)+l)-1

(4 .14)

/eA

The sum over J in the last equation simply corresponds to the number of points in the current filter mask M . If I I denotes the set cardinality, (4 .14) reduces to OF

=2( M -1) .

(Or

(4 .15)

Finally, combining (4 .9), (4 .13) and (4 .15) yields the following result : a3`= MI - 1 .

(4 .16)

Or

Vol . 17 . No . I, April 1992



12

P . Salembier / Adaptive rank order based filters

This equation allows a very simple optimization of the rank parameter . Let us illustrate a possible application of this algorithm . In this example, a constant signal is corrupted by a zero mean noise which has an asymmetrical probability density function . An example of such a signal is given in Fig . 9 . The noise average is null, but the positive deviations are higher and less dense than negative ones . Suppose that a ROF is used to estimate the original constant value . The filter mask is a fixed window of length 15, and the only parameter to optimize is the filter rank . The associated convergence parameter R is set to 10 -7 . Figure 10 presents the evolutions of the normalized rank and the mean square error of the estimation . The initial filter is an erosion and it converges to a ROF with a normalized rank around 0 .7 . As a reference, the level r=0 .5 corresponding to the median filter is drawn . As mentioned in Section 2 .1, this example illustrates a possible use of general ROF . Of course, if necessary (4 .11) and (4 .16) can be jointly used in order to optimize at the same time the rank parameter and the filter mask . Note that in this case, the mutual dependence of the filter mask and of the rank parameter is neglected . In practice, this simplification is not a source of difficulty if the convergence parameters (respectively the time evolutions) are not too high (respectively rapid) . Let us see now how this approach Signal 200-

too-

Noise

k

t'

0

I

Mean value

1

-100

. 0

40

60

80

Time Fig. 9 . Asymmetrical noise. signal eroaroahig

I

Median filter

MSE

0

ltcations Fig. 10 . Rank optimization with a ROF .

can be applied to more complex filters such as combinations of two dual ROES .

5 . Adaptation of more complex ROBFs 5 .1 . Adaptation process

The goal of this section is to describe how the approach used to optimize a single ROF has to be modified in order to deal with combinations of the two dual operators E,, M [ ] and D,, M [ ] . Let us start with a simple combination : suppose that the filter D,, M[E,, M[ ]] is used to produce an estimate yt of a desired process d,, and that r or M have to be optimized in order to minimize the MSE or the MAE . Let us denote by x i the filter input signal . As previously, the explicit formulation giving the filter output as a function of its input is not desirable . Instead, two implicit formulations involving sign functions are used . The first one describes the processing achieved by the first filter E, .,,,[ ] and the second one deals with the dual filter D,, M [ ] . Denote by x ; the output of the first filter :

100

x ;=E, .M[xr1, yi = D,.M[Er.M[ri ]] = Dr .M[x']

(5 .1)



P . Salembier / Adaptive rank order based filters

The implicit function corresponding to the first filter is the same as the one given by (3 .8) . The only modification is to rename the output y, by x ; : F(m;, x,_;, x;, r)

13

derivatives of the two implicit functions F and G given by (5 .3) and (5 .4) . Taking into account that ;/0m,=0 for all j#k, ar/amk =0 and 8x, +,, ;/ am am, =0, the differentiation of G (eq . (5 .3)) gives the following result :

_Y z(sgn(m,)+I) x(sgn(x, ,-x ;)+2r- I)+1-2r =0 .

dG

OG

dm k

am k

(5 .2)

The differences between D,,,,,[ ] and Er,,[ ] are due to the rank which becomes 1-r instead of r and to the reflected mask . This reflection can be written in a very simple way by noting that if a particular location j belongs to the mask M, then -j belongs to the reflected mask . The implicit function describing the filter D,, M[ ] is then the following : G(m,, x ;+ ;,

, ,A

OG

ax' + ;-

ax' + ,-

am,

+G a ay, =0 . a v , am k

(5 .5)

Similarly, the differentiation of F (eq . (5 .4)) yields dF

8F +

dm k

am k

aF ax; +f

ax

;+t

(5 .6)

0

amk

Combining (5 .5) and (5 .6) gives the final partial derivative :

y,, r) +

OG

_ >, z(sgn(m;)+l)

a y;

x (sgn(x ;+ ;--y,) + 1 - 2r) + 2r-1

am,

amk

OG J FA

ax;+,

-(OF / OF am k/ ax ;+

(5 .7)

aG ay.

=0 .

(5 .3)

The last formula involves the output of the first filter E,_,,,[ ] at location i+j' . Substituting x ;+,- for x', in (5 .2) gives F(m,, .x,+;•

r, x ;,,,, r)

It is now necessary to compute all the partial derivatives involved in (5 .7) . Following the same approach as in Section 4 .2, we have successively : From (5 .3), OG S,,.,k(sgn(xl+k-yi)+1-2r),

_ Y, }(sgn(m,)+1)

amk=

( 5-8)

E,7

x (sgn(x, +, -,-x ;+,-)+2r-1)+ 1-2r =0 .

(5 .4)

Having these two implicit formulations of the filtering process, the adaptation problem can be investigated . 5.2. Optimization of the filter mask for Dr,M[Er, .H[ ]]

As in Sections 4 .2 and 4 .3, the partial derivative of the filter output y ; with respect to a particular continuous value Mk associated with the search area can be obtained by computing the total

aG

sgn(m,)+l S(x, ., 2 ax;" - =2 8G ay i

And from (5 .4), OF am,

=0o ;) S(( .)sgnx,+, - k -x+,' +2r-1),

(5 .11)

(OF

(5 .12)

ax , +,Vol . 27. No. 1,

April

1992



14

P. Salembier / Adaptive rank order haled filters

Combining the last six equations gives aY, - , z St

am k

, sgn(xi +

Let us illustrate the properties of this algorithm by several examples .

y ( )+1-2r

sgn(mj ) + 1 j'FA

x (sgn(x,+y'-k-xiy')+2r-I)) . (5 .13) . As Let jo denote the index j' such that =y, this particular location jo has produced the filter output y,, it is certain that sgn(m,,) = 1 . Then, the sum over j' in (5 .13) can be simplified :

EXAMPLES The experiments done in Section 4 .2 for ROFs are reproduced here with the more complex filter D, .w[E, M [ ]] . The first test deals with the filter convergence in the case of a noisy sine wave . Note that it is difficult to talk about a real prediction problem . Indeed, if a filter mask is causal, its symmetric is anticausal . In fact, the computation of y ; always Error

E sgn(m,.)+1 2

St n

1't A

Ii = 0.0001

x (sgn(x;+j'-k-x,'+i')+2r-1)

=sgn(x,+j,-k-x ;+j,)+2r-1 .

(5 .14) 0

Finally, the gradient estimate is a y; am,

µ=

0.001

N=

0.01

z(sgn(x ;+k - yr) 0

+sgn(x,+f,-k-x;_;,)),

(5 .15)

with y,=D, .,M[E,,M[x,]], x'=E,_,,[x,] and jo such that xi+ l,= y, . g As in Section 4 .2, U,,,k has been omitted, and in practice the parameters Mk are adapted between -I and 1 . Note that the rank parameter does not appear in (5 .15) . This means that the adaptation algorithm is the same for an opening, a closing or a two-stage median filter for example . The adaptation formula involves two sign functions, each one corresponding to one of the two filters E,,,,[ ] and D, ,,y[ ] . The term sgn(x';+k-y,) takes into account the error which is introduced by D,,,,[ 1, whereas sgn(x,+1,-k-x;+jo) concerns the error produced by the first filter Er .,,,[ ]- Thus, it seems natural that this second correcting term involves the only point x, ., which is selected by D,,,,[ ] to produce y, (remember that the output y, is always one of the input elements) .

~--~ "

0

2000 Iterations

Fig. IL Mask convergence with a two-stage median filter .

Error

Ii = 0.0001

A = 0 .001

f

- - --r -

. -v -

1

11= 0.01

2000 Iterations Fig. 12 . Mask convergence with an opening .



P. Salembier / Adaptive rank order hosed filters Error

15

Signal

I

p=0.0001 I

11=0 .001

fJh-vIP, IIf .\UfJI IIf A' N ,' V

Closing

Two stage median filler

Opening

I

1 0

p = 0 .01

I \

2000

100

Iterations

Time

Fig. 13 . Mask convergence with a closing .

Fig. 14 . MSE adaptation .

involves the value x, in the first filtering process . This can be seen in (5 .4) for example . The set of original points used to compute xi, and thus y,, is It obviously contains x ; . This remark has an important implication in the case where the desired signal is the original signal itself, d;=x, . In such a situation, there is always a perfect solution for the mask optimization consisting of a single arbitrary point j,, . Indeed, both E,,,,,[ ] and D,,,,[ ] become simple translations and their combination gives the identity

Signal S

1' vV

(5 .16) D, . .M[xi] =y+ =x ; .;p = x,

V

V

s J

\

V

\

V

r

Two stage median filter

\

Opening V

Y

~Y

n

I

Original V

E,, u[xi ] = x,' = z,-io

and

Original

lI

`I V

0

100 Time Fig . IS . MAE adaptation .

Fortunately, it is easy to avoid this trivial solution by constraining the number of points of the solution . For example, if the set {m;} of real values of the search area is such that only one value is above zero, one can shift the whole set of values in order to get at least two positive values . This constraint is imposed in the following tests . Figures 11, 12 and 13 present the evolution of the MSE for a two-stage median filter, an opening and a closing, respectively . Here also, by steps convergence of the MSE is observed .

Finally, Figs . 14 and 15 compare the restoration performance in the case of MAE and MSE optimization, and for various filters . Here also, it can be concluded that the MAE criterion is more robust than the MSE . It is interesting to note that the opening allows for restoration of the negative parts of the sine wave, whereas the closing has a similar effect on positive parts . Of course, two-stage median acts on both parts . Volt 2i, No . I, April 1992



16

P. Salembier / Adaptive rank order based filters

5 .3 . Optimization of the rank parameter for

5 .4 . Optimization of the filter mask for general ROBFs

Dr. .M[Er,M[ ]]

The approach used previously can be applied to optimize the rank parameter . Unfortunately, it gives an estimate of the gradient which is zero . This does not mean that the rank value has no influence on the output, but that the LMS approximation is not suitable for optimizing this parameter . Nevertheless, formula (4 .16) derived for single ROFs was applied in order to evaluate its usefulness in the case of D,,,,,[E,,,,,[ ]] . As in Section 4 .3, a constant signal is corrupted by asymmetrical noise and the rank parameter is optimized to reduce the MSE. The window length is equal to 15 and the convergence parameter A to 10' . Results of this test are given in Fig . 16 . One can see that the adaptation is successfully achieved in a way similar to that of a ROF . The use of this formula can be intuitively justified by the fact that for both filters E,, M [ ] and D,, AM [E, M [ ]], the rank parameter controls a statistical bias in the estimation process (see Properties 1 and 2 of Section 2 .1) and that (4 .16) allows reduction of this bias .

The gradient estimate of the filter output given by (5 .15), involving two sign functions for the composition of two filters, illustrates a general result which can be extended to more complex filters. For an arbitrary composition of ROFs, the approach used in Sections 4 .2 and 5 .2 leads to the definition of a set of implicit functions, each one describing the filtering process of an elementary ROF in the composition . Then, by computing the total derivatives of the functions and by solving the resulting set of equations, the gradient estimate of the filter output can be expressed as a function of several partial derivatives . Finally, the following general rule is obtained, For a particular mask location k, the gradient of the filter output y i with respect to the parameter m,, is proportional to a sum of sign functions plus a correcting term depending on the rank r . Each sign function is dedicated to a specific rank order operator E,,,,,[ ] or D,, .M [ ] . It compares its output, which will be selected by the next filtering process in order to produce the final output y i , with its input at location plus or minus k depending on whether the rank order operator uses the original mask (E,,,M[ ]) or its reflection (D,,m[ ]) . The correcting term takes into account the rank parameter . Its value is equal to 2r-1 times the number of filter E,, M[ ] involved in the composition, plus 1-2r times the number of filter

MSE

D,wJ ] .

Let us describe this rule more formally . Suppose that a general ROBF f,,,M O is constructed by composition of m elementary ROFs, and that thej th ROF in the composition is denoted by a

2000 tteaations Fig . 16 . Rank optimization with a ROBF .

Signal processing

J r,M'( ),

(5 .17)

where r is the rank parameter, M the filter mask and the upper index indicates both the position j



P. Salembier / Adaptive rank order based filters

of the filter in the composition and its type defined by e;=1 if the ROF is E,M[ ] and e,=-1 if the ROF is Dr .t[ ] . The global filter is therefore Yt=.`r.n+(x, )

- J , .,w

rEM 1,( . . . f~M

(I :,

(. . . I

pM (xi))))

17

Let us denote by i+A; the location at level j which corresponds to the output value y ; : xi{ A ,=y ;

Vj, 0

j < m with 4,,,=0 .

(5 .20)

With these notations, it can be shown that the gradient estimate of y, is given by the following formula :

(5 .18) Let us denote by x;" the intermediate signal which is the output of the jth filter :

87th

+'-e". ( s gn (-x;--=,k-

)

+e,(2r-1)) . xv'=.f"b9'`d( . . . f r.

(5 .19)

N

with the two following conventions : Y~

X

'.

X

1

= Yi

As stated previously, it is necessary to follow the trajectory of the final output y; in the various filtering levels . Indeed, the output of a ROF is always one of its input values at some location . Thus, the final filter output at location i denoted by y;=x , 0 comes from a particular location i+A,,, l at the level m-1, that is the value x,"_' A ,", . In turn, this value comes from a particular location i+A,„_ 2 at the level m-2, etc . This procedure is illustrated in Fig . 17 .

Level m

I I I I I I I I I

I

(5 .21)

This rule can be used to easily derive an adaptive algorithm for an arbitrary succession of rank order operators . As an example, the adaptation formula for the filter E,_,u[D,,,,,[D,,,,,[E,,,,,[ ]]]] involving four rank order operators is given by ay` =1

8mk

[sgn(x ;31k-y,)+sgn(x ;+'A,+r, + sgn( z " A,-k - x;+) n,) +sgn(x ;_A,-k-x,+'A,)]-

(5 .22)

This formula is valid, in particular, for open- close and close-open filters . Let us apply now these adaptive algorithms to practical image processing problems .

I I I I

6 . Application to image processing Leveeeel m-1

I I I I I I I I 1 1 1 1 1 1

71' '

ltAm-1 Level m-2 A

I I I I I

I I I I I I I I I

442'

i+A

I _4+4+4+4+1

:1

;

[ (1)

Level 1

t

Level 0 I AO Fig . 17 . Trajectory of the output in the various levels .

The goal of this section is to describe some practical examples which demonstrate the ability of adaptive rank order based filters to cope with textures, noise and nonstationarities observed in images . The theoretical approach presented in the previous sections can be directly used for images when a specific scanning path has been chosen . From the practical point of view, the images in this section are scanned following a zig zag path . This choice is made to preserve a certain continuity or stationarity in the signal and therefore to allow an easier convergence of the algorithm . Moreover, the image is scanned twice. The continuous values associated with the search area are set to arbitrary Vol . E7 . No . I,

April 1992

18

P. Salembier / Adaptive rank order based f lters

values (equal to one, so that the initial filter mask is the search area itself) . Then, the image is scanned following a zig zag path from top to bottom . This first scanning should yield the convergence of the algorithm . Finally, the scanning is inverted (from bottom to top) in order to create the output image with the optimum mask . Note that the filter is still optimized during the second scanning in order to deal with possible signal nonstationarities . 6 .1 . Noise cancellation with adaptive ROBFs

The first example deals with a classical problem of impulse noise suppression . The original image represents a close view of a straw screen [1] and is shown in the upper left part of Fig . 18 . It consists of vertical structures, more or less periodic . The original image is corrupted by impulse additive noise (so-called `salt and pepper' noise) as presented in the upper right part of Fig. 18 . If a

nonadaptive median filter is used to cancel the noise, it will remove both the noise and the vertical thin structures . However, efficient noise suppression can be achieved with an adaptive median filter used as a predictor . It consists in estimating the value of each current pixel with the help of its neighboring pixels . The parameter to optimize is then the filter mask . From a practical point of view, the search area is a 21 x 21 square window excluding the central point which has to be estimated, and the adaptation is performed according to the MAE criterion . As discussed previously, the image is scanned from top to bottom following a zigzag path, and when it has been entirely processed, the scanning is inverted (from bottom to top) . Following this procedure, the filter mask converges towards a vertical structure consisting of ten points on both sides of the point to estimate . The predicted image is presented in the lower part of Fig . 18 . One can see the very good noise suppression

Fig. 18 . Noise cancellation on straw screen . Signal Processing



P . Salembier 11 Adaptive rank order based filters

capability while preserving the background structure . The following two examples show the system ability to deal with nonstationary signals . Note that there is no benefit in scanning the images twice . Thus, the results presented here have been obtained after a single scanning . The original images are presented in the upper left part of Figs . 19 and 20, respectively. The first one represents a close view of a tree stump [1] with a wood knot, whereas the second one shows a man speaking in front of a curtain with vertical structures . These original images are corrupted by additive impulse noise as shown in the upper right part of Figs . 19 and 20 . In the case of Fig . 19, an adaptive twostage median filter is used in order to restore the image . The search area is a square of size 11 x 11, and the initial mask is set to a horizontal line of length equal to 11 . The image in the lower left part of Fig . 19 shows the restored image when the

19

convergence parameter is very low, A.=10 - s (MAE adaptation) . As can be seen, the parts where the wood structure is horizontal are well restored . But, all details have been removed from the area corresponding to the wood knot . Indeed, the convergence speed being quite slow, the filter mask is not sufficiently modified inside this restricted area . In contrast, if the convergence speed is higher, the algorithm can actually match the filter mask to the new structure . This is illustrated by the lower right part of Fig . 19, which is obtained with a convergence parameter A of 10 -3 . Similar conclusions are obtained with the image of Fig . 20 . Here, noise cancellation is achieved with adaptive median filter . The search area is a square of size 7 x 7 and the convergence parameter is set to 10 ` (lower left part) and 10 - ' (lower right part), In the case of low convergence, most of the noise has disappeared but some details have also been removed . Note in particular the spectacles of the man or his hands .

Fig. 19 . Noise cancellation with nonstationary signals . Vol . 27, No 1 . .April 1992

20

P. Salembier / Adaptive rank order hased filters

Fig . 20 . Noise cancellation with nonstationary signals .

In fact the filter mask has taken a vertical shape corresponding to the background and the system is not able to rapidly modify the mask shape for other small regions . On the contrary, if a larger value of the convergence parameter is used, the filter is possible to react rapidly in case of nonstationarities and both remove noise and preserve most of the details . 6 .2 . Image restoration with adaptive ROBFs Adaptive ROBFs are attractive solutions to perform certain types of image restoration . Suppose for example that an image has been corrupted in such a way that most of the original image is known but some samples are missing . Furthermore the locations of the missing samples may or may not be known . An example of such a corruption can be found in image communication where transmission errors may create missing samples . Signal Promssing

Another example is the case where the image has been damaged by dirts or by human actions . The left part of Fig . 21 shows two images which have been corrupted by black lines . In fact this problem is quite similar to the noise cancellation problem which has been illustrated above . The major difference is in the random nature of the corrupting process . In the examples of Fig . 21, the missing samples are not randomly located even if their precise positions are not known . Let us describe the restoration of these two images . In both cases, we know that the missing samples appear as black pixels . Therefore, the restoration can be performed by a closing which preserves bright components and removes black ones . The perturbation can be viewed as an impulse event, thus the MAE criterion is used to optimize the filter mask . For the first image, the search area is a square of size 11 x 11, the convergence parameter is set to 3 x 10-4 . In order to be sure to remove the



21

P . Salembier / Adaptive rank order based filters

Fig.

21 .

Image restoration with adaptive ROBFs .

black lines, the mask is constrained to have at least 18 points . As mentioned in Section 5 .2 this constraint can be very easily imposed by shifting the whole set of continuous parameters of the search area . The constraint is necessary because the image is made of rather smooth but non-constant surfaces . This means that a very good prediction can be achieved with a very small number of points . If the mask contains only a few points, the filter may not perfectly remove the black lines . For comparison purposes, let us mention that the lines width is roughly equal to 4 pixels . This means that if a fixed closing is used, the size of its mask should be larger than 4 x 4= 16 points . The upper right part of Fig . 21 shows the restored image . As can be seen, the quality is quite good . It is difficult to locate the presence of the perturbation . Only some edge points are not perfectly restored. In this example, the mask of the filter has converged toward a more or less isotropic shape reflecting the absence of

specific orientation in the original image . The next example illustrates the case where the original image is textured . The image presented in the lower left part of Fig . 21 consists of a canvas view which has been corrupted by two black lines . The restoration process should remove those lines, but it should also restore the background texture . For this application, adaptive techniques are very attractive since they allow extraction of the main texture characteristics . As previously, the restoration is achieved with closing optimized following the MAE criterion . The search area is rather large (size= I I x 11) and the filter mask is constrained to have at least 20 points . The restoration result is shown in the lower right part of Fig . 21 . As can be seen, the black lines have been successfully removed and replaced by the `original' texture . Note that the restoration process has also partially removed the noise and the defect which are present on the Vol . 27 . No I

April 1992

22

P. Salembier / Adaptive rank order based filters

original image . A possible solution to avoid this phenomenon is to process the image only for regions which are close to the black lines . 6 .3. Detection of texture defects This last experiment deals with the smoothing effect of adaptive two stage median filter (Dr,M[E,,M[ ]], with r=0 .5), in the presence of textures . Three original images are presented on

the left side of Fig . 22 . They represent views of a piece of wood, of a wire mesh and of a piece of textile, respectively . In each case, the background texture exhibits some anomalies . For the piece of wood, the main defect is caused by an impact on the surface and there are some other small dark points due to dust . Regarding the wire mesh, the texture anomaly corresponds to a deformation of the mesh by a circular tool . Finally, the textile surface exhibits a weaving defect .

Fig . 22 . Texture defect detection . Signal Processing



P. Salembier / Adaptive rank order based filters

The mask of an adaptive two-stage median filter is optimized with the aim of smoothing the images . In all cases, the search area is an entire square of size 13 x 13 . As our intention is to have a noticeable smoothing effect, a lower bound on the number of points of the filter mask is defined . Its value is equal to 15 in this test . This constraint is necessary, since the optimum unconstrained filter mask may comprise much less points (having of course removed the trivial solution involving an arbitrary single point as discussed in Section 5 .2) . The same scanning procedure as in the previous test is adopted . The results are presented in the middle column of Fig . 22 . It can be seen that the smoothing effect allows noise reduction with a very good preservation of the background texture . On the specific texture anomalies, the filtering effect depends on the background . In the case of wood surface, the defects are removed and the filter performs a vertical interpolation, whereas the hole in the wire

23

mesh disappears and is replaced by the `original' mesh structure. Finally, the weaving defects have been partially removed and replaced by the hexagonal periodic texture . In the case of wire mesh and canvas images, the filter mask has converged toward a disconnected shape which reproduces the background periodicity . Adaptive prediction can be used as feature extraction for texture defect detection . Indeed, a defect can very often be considered as an impulse and unpredictable event in the image . Therefore, adaptive ROBFs are an attractive solution to extract the predictable part of the original image . A low prediction error characterizes a defect-fret area, whereas a high prediction error in absolute value reflects the presence of a defect . The images of absolute prediction error are shown on the right side of Fig . 22 . As can be seen, the defects are enhanced and appear as bright peaks in rather low level noise .

Fig . 23 . Smoothing effect of various adaptive ROBFs . vol . 77. No . I . April 1992



24

P. Salembier / Adaptive rank order based filters

The last example illustrates the difference of smoothing effects of three adaptive filters, namely two-stage median filter, opening and closing . The original image is the canvas image presented previously . The upper right part of Fig. 23 gives the result obtained with a two-stage median filter . The search area is a square of size 13 x 13 . As in the previous examples, it suppresses noise and restores the `original' visual appearance of the texture . The opening (lower left part of Fig . 23) produces a `lower estimate' of the original image . It has an effect on only those image pixels with brightness exceeding the predicted value, and the dark parts of the defect remain unchanged . As expected, the closing (lower right part of Fig . 23) has a dual effect .

7 . Conclusions In this paper, the problem of adaptive optimization of a class of nonlinear filters referred to as rank order based filters has been addressed . This class includes, as special cases, median, classical rank order and morphological (with flat structuring element) filters . Its definition highlights a common origin based on sorting operations . The parameters which can be optimized are the rank order involved in the sorting operation and the filter mask . The optimization intends to minimize the mean square or the mean average error between the filter output and a desired signal . It is achieved by a recursive adaptation of the current solution similar to the least mean square algorithm . The optimization problem comes up against two main difficulties . The first one deals with the sorting operation(s) involved in the filtering process . Indeed, the derivative of such an operation is difficult to obtain or manipulate . Therefore, an alternative equivalent formulation of rank order filtering has been proposed . Its main advantage is that it involves only sign functions which are easily derivable . This expression is obtained using an implicit formulation of the filtering process . Signal Processing

The second problem concerns the binary nature of the filter mask which makes it difficult to adapt by a recursive approach . To solve this difficulty, a search area which should contain the optimum filter mask has been defined . A continuous parameter is associated to each possible location of the search area, and the filter mask is defined by thresholding this set of continuous values . The optimization problem is then redefined as the optimization of the set of continuous parameters instead of the direct optimization of the filter mask . This approach leads to very simple algorithms which can be used for any combination of rank order filters . The convergence properties and the difference between MAE and MSE optimization have been illustrated on some practical examples with one-dimensional signals . Finally, applications to image processing were presented . They show the ability of adaptive nonlinear filtering in coping with different kind of textures, in removing noise and in dealing with nonstationary signals .

References [1] P . Brodatz, Textures - A Photographic Album for Artists and Designers, Dover, New York, 1966. [2] E .R. Dougherty, "Minimal search for the optimal meansquare digital gray-scale morphological filter", Proc. SPIE Visual Communications and Image Processing 90, Vol . 1360, Lausanne, Switzerland, pp. 214-226 . [31 M . Gabbouj and E .J . Coyle, "Minimum mean absolute

error stack filtering with structural constraints and goals", IEEE Trans . Accost . Speech Signal Process ., Vol . ASSP 38, No . 6, 1990, pp . 955-968 . [41 J.H . Lin and E .J. Coyle, "Minimum mean absolute error

estimation over the class of generalized stack filters", IEEE Trans . Acoust . Speech Signal Process ., Vol . ASSP 38, No . 4, 1990, pp . 663-678 . [5] J .H . Lin, T .M . Sellke and E .J . Coyle, "Adaptive stack filtering under the mean absolute error criterion", IEEE Trans. Acoust . Speech Signal Process ., Vol . ASSP 38, No . 6, 1990, pp . 938-954 . [6] P . Maragos and R.W. Schafer, "Morphological filters Part I : Their set-theoretic analysis and relations to linear shift-invariant filters", IEEE Trans . Aeoust_ Speech Signal Process ., Vol . ASSP 35, No. 8, 1987, pp . 1153-1169. [71 P . Maragos and R .W. Schafer, "Morphological filtersPart 11 : Their relations to median, order-statistic, and stack filters", IEEE Trans . Acoust. Speech Signal Process., Vol . ASSP 35 . No . 8 . 1987, pp . 1170-1184 .



P . Salembier / Adaptive rank order based filters [81 T.A . Nodes and N .C . Gallagher, "Median filters : Some modifications and their properties", IEEE Trans. Acoust . Speech Signal Process . . Vol . ASSP 30, No . 5, 1982, pp. 739 746 . [9] 1 . Pitas and S . Vougioukas, "Adaptive nonlinear filters based on order statistics", Pear . EUSIPCO 90, Barcelona, Spain, pp . 397-400 . [10] P . Salembier . "Multiresolution decomposition and adaptive filtering with rank order based filter -Application to surface defect detection", Proc . Internal . Coal. Acoust. Speech Signal Process ., 1991, Toronto, Canada . pp . 2389-2392 . 111] P. Salembier and M . Kunt, "Size-sensitive multiresolution decomposition with rank order based filter", IEEE Trans . Pattern Anal . Machine Intell . . submitted .

25

[12] K . Saarinen and Y . Neuvo, "Adaptation of weighted median filters", Proc. EUSIPCO 90 . Barcelona, Spain, pp . 269-272 . [13] J . Serra, Image Analvris and Mathematical Morphology, Academic Press, New York, 1982 . [14] P .D. Wendt, E .J . Coyle and N .C . Gallagher . "Stack filters", IEEE Trans . Acoust . Speech Signal Process ., Vol . ASSP 34, No . 4, 1986, pp. 898-911 . [15] B . Widrow and S .D . Steams, Adaptive Signal Processing, Prentice Hall, Englewood Cliffs, NJ . [16] L . Yin, J . Astola and Y . Neuvo, "Adaptive stack filtering with application to image processing" . IEEE Trans . Acoust . Speech Signal Process ., submitted .

voi . 27 . No . t . April 1992