The effect of median filtering on edge location estimation

COMPUTER GRAPHICS AND IMAGE PROCESSING 15,224-245

(1981)

The Effect of Median Filtering on Edge Location Estimation* G. J.

YANG+

AND

T. S.

HUANG$

School of Electrical Engineering, Purdue Unioersi@, West Lufwette,

Indiana

47907

Received March 25, 1980 We study the effect of noise reduction preprocessing, specifically median filtering and averaging, on the accuracy of edge location estimation using least squares. The original edge is either a step or a linear ramp, corrupted by white Gaussian noise or binary symmetrical channel noise. The surprising conclusion is that in the case of white Gaussian noise, neither median filtering nor averaging improves the estimation accuracy. In the case of binary symmetrical channel noise, median filtering does improve the estimation accuracy for ramp edges, which are reasonable models for real-life edges. 1. INTRODUCTION

An important issue in image analysis by machines is the effect of preprocessing (such as noise reduction, contrast enhancement, and restoration) on further machine analysis (mensuration, detection, and classification) of the image. Many preprocessing methods may improve the subjective quality of images but m fact degrade the performance of further machine analysis. In this paper we present results on a particular aspect of this important issue: The effect of median filtering on the accuracy of edge location estimation using least squares. Some results on the effect of median filtering on shape analysis are presented elsewhere [I]. Compared to the long history of linear filtering, median filtering is relatively new. In the early seventies, Tukey [2,3] was among the first to suggest the use of median filters for signal smoothing. The running median filter was found to be able to keep signal jumps in the signal while reducing noise. More recently, Rabiner ei” a/. [4] and Jayant [5] applied median filters to speech processing. Pratt [6] used median filters to suppress impulse noise in images. Frieden [7] used median filters to reduce the spurious oscillations around the edges after a linear restoration technique was applied to one-dimensional images. Narendra [8] and Shamos [S] approximated the two-dimensional median by the median of the medians in each row of the window to facilitate median circuit design for 2-D image processing. Eversole, et al. [lo], designed IC circuits for median filtering using VLSI technol-. ogy. Huang, et al. [ll] developed a fast algorithm for two-dimensional median filtering, which is based on histogram modification. Although median filters are useful and getting popular, theoretical results on their behavior are scarce in the open literature. Recently, Tyan [ 121 obtained some analytical results for the fixed points of the running median and Justusson [ 139 studied the statistical behavior of median filters. Velleman [14] characterized them by their responses to zero-phase single-frequency sine waves. Motivated by the fact *This work was supported by the U.S. Army Research Office under Contract DAAG 29-79-C-0200 (Technical monitor: Dr. William Sander, ARO; Scientific Liaison Representative: J. J. Fasano, USA ARRADCOM). +Now at Xerox, El Segundo, California. *Now at Coordinated Science Laboratory, University of Illinois, Urbana, Ill. 61801. 0146-664X/81/030224-22$02.00/0 Copyright 0 1981 by Academic Press, Inc. All rights of reproduction in any form reserved.

224

MEDIAN

FILTERING

-

noise linear Step

no

blurring

filtering

Least

or

edge

or

00...01..1

no

-

blurring

.---

square

averaging olmedian

-

225

AND EDGE LOCATION

estimation

filtering

Edge -r

location estimate

~

1

FIG. 1. Blockdiagramfor edge location estimation.

that median filtering tends to preserve edges in images, we investigate quantitatively the effect of median filtering on one-dimensional edge location estimation. Edge location estimation is an important part of the image mensuration problem because edges carry object boundary information. Burnett and Huang [ 151 described a MAP (maximum a posteriori probability) estimate of one-dimensional edge location. In the model used, an ideal step edge is blurred by a linear shift-invariant system and degraded by a noise source. We adopt the same model but search for the least-squares estimate (instead of MAP) of the edge location. Usually, the edge transition probability is unknown, and the MAP estimate is reduced to a maximum likelihood estimate. For the white Gaussian noise case, the maximum likelihood estimate is the same as the least-squares estimate. Without loss of generality, we assume the two levels of the step edge to be 0 and 1. The block diagram of our procedure is shown in Fig. 1. We consider two cases: Case I -no blurring, Case II-linear blurring. For each case, we studied two different types of noise: white Gaussian noise, and noise due to transmission of the signal by PCM through a binary symmetrical channel (BSC). Theoretical results for Gaussian noise are presented in Sections 2-8. Supporting experimental results are in Section 9. For BSC noise, no theoretical analysis was done. Only the experimental results are presented in Section 10. 2. CASE I WITHOUT

FILTERING

In Case I, the observed sequence is the step edge plus noise. We perform no filtering, averaging, or median filtering before estimation to compare the effects of these operations as far as the edge location estimation is concerned. For no filtering, we can model this estimation as a Mary detection problem in statistical communication theory. The set of M vector signals of dimension M is s,= (l,l)...)

l),

s, = (0, 1, . . . ) l),

(1) s,=

(0,O )...) 1).

They are linearly independent and separated by Euclidean distance of

d(Si,Sj) = qT=J

.

For white Guassian noise, the least-squares estimation procedure is equivalent to the maximum likelihood decision rule, which gives the minimum probability of

226

YANG

AND HUANG

error, p,. In the decision space, the M decision regions can be expressed as follows: Ri=

;

i= 1,2 )..., M,

R,,

1’2)

j-1

i#i

where R.. is the left half space (including A’;.) whose boundary bisects the hne L segment SiSj Let the correct decision be Sk. The probability of deciding Sk+,, is P(Sk,nISk)

= 1

P(ZISkW Rk+n

where Z is the observation vector. Due to symmetry, PCS,-,I$)

= P(S,+,IS,).

The total probability of error is P, =

c

P(Sk+nISk)

I-Ksn
The probability of correct estimation is

P, = P(SkISk). For A4 > 2, it is in general not possible to express P, in terms of the composite error function, Q. For a simple case where M = 3, K = 2, n = + 1, we can draw S,, S,, and S, in the plane of xs = 1, as shown in Fig. 2.

= Q’(-+-).

S2(0,1,1) S3(0,0,1) I+=X2

R2

R3

RI

.

S,(l.l,l,

t Xl

Ro. 2. A simple decision space for edge location estimation.

MEDIAN

FILTERING

AND EDGE LOCATION

227

For step size h, PC’ Qf-g-) P, =

P(S,IS,)

=

WS,IS2)

= (1 - e(&-))I +

= 1 - (1-

P(%IS2)

Q(k))?

(7)

When A4 2 4, an exact expression is not available, but some bounds can be derived. For any k, k # 1 or M (meaning that edge jumps do not occur at end points),

[I-

< Pc<[ 1 - Q(%,]'.

Q(&)l"-'

(8)

The reasoning is as follows. For a fixed k and any n such that 1 I(k+nork-n)

IM,

we have (9) Let n be positive in Fig. 3. Define v,j = Si - 3, 1 I i,j I M. Then V k,k+n

to,.

,O,l,

1=

co, *

t t kth (k + n)th column column 1 1 ,O,l)..., l,l,O ,..., 0)

1=

co,.

. . . . . . . . , O,l,O ,...) 0)

=

Vk,k+n+ V k+n,k+n+

s

. ..) l,o )...) 0)

k+n

FIG. 3. Decision space for computing

k+n,k+n+l

error bound.

(10)

228

YANG AND HUANG

Define +Jk+n,k+n+l

=

L(Vk,k+n,

= tan-‘-

‘k,k+n+l)

1

n : natural numbers.

v5 Hence 0” < + k+,, k+n+, I 45”. Since r+ < 90”, the (n + 1)st bisecting hyperplane, which is part of decision boundaries, retains more than (1 - Q(h/2a)) x 100% of the probability computed when only up to n bisecting procedures are done. The same argument holds if we replace Sk +n by Sk_ n and Sk+,,+ , by Sk-,, _, . The total number of bisections is M - 1, and therefore the lower bound of PC is [I -Q(h/2a)l”-‘. The upper bound is [I - Q(h/2a)12 because the bisecting due to with Sk give exactly that probability. In summary, for M 2 3 and sk-l and sk+l assuming the jump does not occur at either end, we obtain

[l - Q($)l”-’

< PC+

- g(%jll.

(I U)

In looking into the distribution of errors, we observe that the probability of error is mostly contributed by small ] n ] because Q(m h/28) decreases rapidly as n increases. It is interesting to note that (i) PC decreasesas M increases, suggesting that the total number of data points should be as small as possible (and still contain the edge transition points). (ii) While the derived bounds are rather loose for large M, they nevertheless suggest that the probability of correct estimation (and hence of error) depends mainly on h/u, the neighboring signal distance in decision space to noise ratio. We therefore concentrate on measuring this signal-to-noise ratio before and after each filtering. 3. CASE I WITH AVERAGING

The averaging process is linear. It blurs the step edge and reduces the noise variance. For an odd window size of W, the resultant edge signal contains a linear ramp with (W - 1) equally spaced points between the two levels. The noise is still Gaussian but correlated. The least-squares estimation no longer gives the maximum likelihood estimate. We do not intend to derive PC for either estimation scheme in this section. However, it is probably approximately true that the resulting PC still depends mainly on &/a, in addition to W, where h, is the Euclidean distance between the blurred version of the two neighboring edge signal vectors, and ua2equals u’/ W. The noise correlation is not considered here, but will be in Section 8. To compute h,, assume the blurring of the edges does not extend to the endpoints. From Fig. 4, which shows two neighboring blurred edges, the

MEDIAN

FILTERING

AND EDGE LOCATION

229

points

FIG. 4. Two neighboring edges.

distance is 2

.w

h

=-

cw

.

(12)

Hence,

(13)

The net result is that the average filtering compare this ratio with that in Section 4. 4. CASE I WITH

does not change this ratio. We shall

MEDIAN

FILTERING

A step edge is a fixed point of the median filter, i.e., it is unchanged by the median filter. Due to nonlinearity of the filter, the median-filtered noisy step edge sequence no longer has the step function as its mean sequence, nor does it have equal variance across the edge. Closed form expressions for the mean and the variance across the edge are not possible because they are expressed in terms of integrals of the error function. I?ecently, Justusson [13] computed such sequences as a function of h/u. As h/u increases, the mean at k points away from the edge approaches a constant S measured in (I, the constant being exactly the expected value of the [(W + 1)/2]th largest order statistic of ((W + 1)/2 + k) independent, identically distributed, unit-variance Gaussian variables, where 1 I k I (W 1)/2. For W = 3, S = E{X,:,} = 0.56 where X2:2 denotes the largest of X, and X2. This is good for h/u 2 3 approximately. Therefore, for window size 3 and h/u 2 3, considering the class of the noisy edge signals (including the ramp edges as will be clear in Case II) and considering the mean sequence of the noisy signal sequence as the input and output of the median filter, this nonlinear filter can be modeled by a linear space-invariant system if we measure the signal amplitude in terms of u. Let the filter coefficients be a,, a2, and a J; then the equations for them are 0.(a, + u2) + ha, = 6, O-a, + h(a, + Us) = h - 6, a,+ a,+

a3= 1.

(14)

YANG AND HUANG

230

The solution is 6 h

a,=a,=-=-

0.56

h/u ’

Considering the mean sequence as the blurring model for median filtering, we can now calculate h,/u,. h, is the distance between the following two vectors, for W= 3 and h/o 2 3:

[O,. , . >O,S,h-

S,h . . . . ..,h]

[O...

S,h -- S.h ,...,

and . ..) 0,

h]

Hence hZ, = 6’ + (h -- 26)’ -t 6’ and ,,,,=,[(I

-Fr+2(;r]“’

+(+$)]“2

zh[l-(~-3$)-2$]

Note that u, is not constant. For the two points around the edge transition,

For other points,

a, =

Therefore ~h %I

_ h[ 1 - 2(8/h) -

+ 6*/h’]

au

1 - l.l2(u/h)

(a = 0.82 or 0.67)

+ 0.31(a2/h2) a

]

($3).

(20)

MEDIAN FILTERING

AND EDGE LOCATION

231

The factor in the brackets 1 - l.l2(o/h)

+ 0.31(+12)

a

> 1 -

if h/u 2 6, where 3.1 I b 5 5.9. Therefore for h/u 2 6 we have (21) Define a figure of merit F for a filter as follows: (24 We have, for median filtering,

F, =

1 - l.l2(u/h)

+ 0.31(u/h)*

(23)

a

As h/u increases, Fm+-,

1

(24)

a

where l/a = 1.22 for the two points around the edge transition, and l/a = 1.49 for all other points. So, asymptotically, there is approximately 20 to 50% increase in signal distance to noise ratio. Recall that for averaging, F, = 1. We summarize the situation below. Median-filtered output no longer has Gaussian-distributed noise (but asymptotically Gaussian as the window size increases). Noise distribution functions and noise correlations are also different across the edge, and the noise sequence is nonstationary. This makes ambiguous the correct value of a,,, and therefore less accurate the performance measure of median filtering for edge estimation. The probability of correct estimation by the least-squares method, PC, also can drop from that of the maximum likelihood estimate because our least-squares estimate does not take noise correlation into account and, due to non-Gaussian distribution, would not be a maximum likelihood estimate even if the noise sequence were truly uncorrelated. 5. CASE II WITHOUT FILTERING

We use averaging of window size Q to model linear blurring in Case II. Without filtering, we have

h no - hi* U no

. U

(25)

6. CASE II WITH AVERAGING

With averaging of window size W, the resultant blurring can be expressed as the convolution of two rectangular weighting functions. The result is a triangular

232

YANG

AND HUANG

weighting function with center coefficient 1/IV and base width 2W points, letting W = Q for simplicity. We can compute the difference of two neighboring edges as follows: s,=

i $2 w2

w2”“‘-9

W

j+...’

2 1 --w2Y w2

+ [step edge at location k] I

- 1z

1

2 W 2 I 9 j+ P.. * 3 F Y*. * 9 --w* 9 w* I * [step edge at location k -t 11

1 2 W ---@9--+‘...‘-+‘...’

2 1 --w2Y w2

* [impulse

at location k with strength h 3

I where

l

denotes convolution.

We then have

ha = Isal

= ---$( W(2W2 + 1)/3)l”-h W= 2,3,...,

=

(26)

d

2 h 30

= 0.8165

6’

ha/% F”=h,,lo,, = 0.816@,

(2’7)

where Q = W = 2,3,. . . Thus, for Q = 3

F, = 1.413. Again, this number is to be compared 7. CASE II WITH

(28)

with F,. MEDL4N

FILTERING

Consider a ramp edge with (Q - 1) transitional points corrupted by white Gaussian noise. Median filtering is applied to it. We are interested in the mean sequence and variance sequence of the output. For window size 3, some of the previous results for Case I can be used. At the points 1 and Q + 1 as shown in Fig. 5, the means are off the ideal value toward the transitional points by a certain portion of h/Q (instead of h).

MEDIAN

FILTERING

AND EDGE LOCATION

233

Q+' Q

FIG. 5. Ideal

rampedge with blurring

width of Q points.

The variances are obtained accordingly. The means for the transitional points are unchanged due to symmetry of the transitional levels in each slicing window. The variances, however, are expected to be larger than those at points 1 and Q + 1. While it is possible to derive a closed form expression for the density functions of the points across the edge after median filtering of any window size W, the result is complicated and to find the means and variances still requires numerical integration. We therefore performed a computer simulation for Q = 3,5, W = 3,5 and h/u = 1,3,5. We plotted two of the mean sequences(based on 400 noisy edges for each plot) in Fig. 6. They agree with the above theoretical prediction. We are now to derive the ratio h, /u,,, for Q = W = 3 and h/o = 5. The ideal ramp edge is replotted in Fig. 7. The mean sequence is 0 ,..., O,& f h, f h,h - &h,.. .,h. The standard deviation sequence is, from the simulation result, . . ,0.67u. Hence 0.670,. . . , 0.67u,0.73u,0.84u,0.73u,0.67u,.

=-

h 1

v3 N-

(32 ‘P

l-4:+12?

h

l-2:+4$,

VT = au %-I

a = 0.67,0.73,

1 1

(29)

or 0.84,

and h

hm -= %I

[ 1 - 2(6/h) a

UV3

F = 1 - 2(6/h) m

+ 4(a2/h2)]

3

+ 4(a2/h2) a

=

1 - l.l2(u/h)

+ 1.25(u2/h2) a

With h/u = 5, F, = 1.231 for a = 0.67, 1.131 for a = 0.73, 0.982 for a = 0.84.

for a 2 3.

(30)

234

YANG AND HUANG 160 eooo

1~¶.000

130 IO00

1ts.000

100 1000

85.0000

70 *oooo

ss.0000

1

90 .oooo

s

10

1s

20

n

30

3s

-40

a 160 so00

1rs1.000

130.000

1ts.000

100~000

05.0000

70 *DO00

ssR.oeoo

b FIG. 6. Mean sequences of noisy ramp edges after median filtering. (a) Q = 3, W = 3, h/o 1 5. @I) Q-5,

W=3,h/a=5.

MEDIAN FILTERING

AND EDGE LOCATION o

o...o-

235

h

0 0 o-o-*.0

0

FIQ. 7. Ideal ramp edge with blurring width of three points.

Asymptotically,

F,+

A = 1.49

for a = 0.67,

1.37 1.19

for a = 0.73, for a = 0.84.

(31)

Compared to the result in Section 6, F, is worse than F, except for the very high signal-to-noise ratio and the optimistic selection of a. Therefore, for the example of W = 3 median filtering performs worse than average filtering for a ramp edge. Note that we had the opposite result for a step edge. Another comment is in order. We mentioned a linear model for median filter of window size 3. It is also valid for the inputs of noisy ramp edges. A direct convolution of an ideal ramp edge with the derived filter coefficients will justify this fact. The remarks in Section 4 regarding median filtering also apply here. We summarize the results for F measurement in Table 1. 8. ERROR PROBABILITIES

In this section, we shall derive an approximate formula for P,. Assuming a step edge as in Section 1 and assuming the true edge jump occurs at the kth point and k 2 n do not exceed the signal interval, we can readily show that

(Sk,” - S,)*(S,-, - S,) = 0 and

d(&,S,,,)

= fi

for all n

TABLE 1 F Measurement for Edge Location Estimation Step edge Averaging, any window W

Ramp edge

1.0

0.816@ (Q = W) 1.413(Q = 3)

1.21 - 1.49

1.19 - 1.49(Q = 3)

Median, W = 3,

high h/a

(32)

236

YANG AND HUANG

‘k-n

5

k+n

FIG. 8. Decision space for computing error.

Thus in the plane of Sk, Sk+,,, and Sk-,, these three points form a right equilateral triangle as shown in Fig. 8. Define dj = d(X, sj) where X is the observed data vector. The probability of error given that Sk is the true edge signal is P(EIS~) = P(any dj < dk,j # k)Sk) = P(d,+,


5 P(4+1

< d,ord,_,<

+ P(d,+, + ... . Let P, = P(d,+,

dk)

< d, or d,-,

< dk)

(33‘)

< d, or dken < dk). From Fig. 8, it can be shown that P,=l-Q2 =24x3

-Gh ( 1 7

-4q.

(341

The sequence of P, decreases rapidly when the signal-to-noise ratio is not too low. For example, if h/u = 3, then P, = 0.129, P2 = 0.033, P3 = 0.009.

For h/u = 4,

PI = 0.045, pz = 0.0048. Therefore, we can use P(dk+l < dk or d,-, < dk) or P(dk+, < dk) + P(d,-, as an approximation of P, for h/u 2 3. Geometrically, we can easily see that dk+, < dk is equivalent to

< dk)

MEDIAN

FILTERING

237

AND EDGE LOCATION

The left-hand side is a linear combination of zero-mean Gaussian variables, and hence also a Gaussian variable with zero mean. Its variance changes with the noise correlation. Equation (35) is useful for the case of averaging of noisy edges. After averaging of window size W, the noise autocorrelation is

,=’

w( = 0,

*A!!

111-c w,

W’1

otherwise.

(36)

The edge component becomes blurred as described before. We are to compute the variance of the LHS (left-hand side) of (35): LHS = N. (-Sk + Sk,,) Wk9Sk,,) = (n,,.*-,n,,-,n,)* 1

(

h/m

0 )...) 0,s

)..., $

,...) -$,o

)..., 0

t

k th column =- -1

5

\r

‘k+iy

i---p

where 2p + 1 = W. Var[LHS]

5

=$E i--p

=+

(

i

wu;+2

+(w-z)q+;+...

2 1: o,‘- w 3 c-0 2

3

where aa = e’/ W.

2

’

nk+ink+j

j--p

[

(W-l)

w- w 1) ff,’ +2&r

+ -Lo2

WB I)

(37)

238

YANG AND HUANG

For the RHS (right-hand side) of (35)

Therefore, a performance measure of the step edge location estimator after averaging is

and a figure of merit for averaging is F, = F, decreases as W increases and is always less than unity. This is more accurate

than our previous calculation, where we did not take noise correlation into account and obtained u,’instead of u,’f Win E@ (37). In the computer experiment, we have M=

100,

h = 150 - 50 = 100,

(I = 31.6. The predicted error bound is, for no filtering on the step edge, P,=2Q

i

;i

1

,

= 2Q(1.58), = 0.1142. From 100 experiments, the observed error probability is P, = 0.07. For averaging with window size 3,

Pe=2Q

--

= 2Q(1.117) = 0.264.

From 100 experiments, we get PC= 0.27.

MEDIAN

FILTERING

I -80

-6:

FIG. 9. The probability filtering.

-4G

-20

239

AND EDGE LOCATION

0

density of VJ for computing

t 2G

I

I 45

I

edge location

60

80

estimation

n error after median

For median filtering, it is again difficult to derive a general and accurate expression. LHS of (35) with k + 1 being replaced by k - 1

Define this random variable as 77.Using a previous result, N*(S,-,

- Sk) = an,-,

+ (h - 26)n,-,

+ f3flk,

where 6 = 0.560 for h/u 2 3. k(N-(Sk-,

- Sk)) = E(b,-,

+ (h - 26)n,-,

+ 6n,)2

(9

The difficulty is in expanding Eq. (40) since ni is a nonstationary and nonGaussian sequence and is correlated in an unknown manner. Also, the distribution of N.(S,- i - Sk) is not Gaussian and therefore the Q function does not give the probability of error. We can, however, perform computer simulation for the numerical solution of a, and the distribution of g. We used 100 noisy edge sequences in simulation, and h = 100, u = 31.6, k = 50. In addition, 6 = 0.56~ was used. The resultant probability density function was plotted in Fig. 9. The standard deviation is 29.75. Also plotted in Fig. 9 was the Gaussian distribution with zero mean and standard deviation 29.75. Using Eq. (16) we have d(s,-,,S,)=

h,=

69*28

240

YANG

AND HUANG TABLE

2

Edge Estimation Pcrformanu! under Gaussinn-Noise (h/a = 3.16) step edge

No filtering Averaging w=3 Median filtering w-3

Ramp adee w=

points)

p.

Errorsizc std. dev.

0.07

0.39

0.39

0.27

0.59

0.49

0.75

0.21

0.60

0.50

0.86

Therefore, predicted probability

Error size std. dev. _.... _ 0.66

P, --..-

of estimation error is, using Fig. 9,

p, = v47l

> kn/2)

= 2P,(Tj > 34.64) = 2 x 0.10 = 0.20. Computer experiments sequences gave

of least-squares

estimation

using the same 100 noisy

P, = 0.21. This is quite close to 0.20, which is the predicted value by Eq. (35). If we regard TJ as Gaussian with zero mean and standard deviation 29.75, then the predicted probability of estimation error is PC= 2Q = 2Q(1.164) = 0.244. This is farther away from 0.21 than the simulated result and agrees with the fact that TJis not Gaussian-distributed. TABLE Distribution Error size

-3 -2 -1 0 1 2 3

3

of Error !Sizein Table I

stepedge

-pedge

No fib.

Aver. fib.

Mad. fililt.

No fililt.

Aver. fib.

Mad. filt.

0 2 4 93 0 1 0

0 0 15 71 12 0 0

1 0 9 79 10 0 1

0 1 14 61 23 1 0

0 1 21 51 25 2 0

1 2 20 50 25 1 1

MEDIAN FILTERING

AND EDGE LOCATION

241

Cd)

(e)

(cl

(f)

FIG. 10. Typical edge sequences.(a) original step edge with Gaussian noise, SNR = 1Odb.(b) 3-point average filtering of (a). (c) 3-point median filtering of (a). (d) original 3-point ramp edge with Gaussian noise, SNR = 1Odb. (e) 3-point average filtering of (d). (f) 3-point median filtering of (e).

242

YANG AND HUANG TABLE 4 Edge Estimation Performance under Binary Symmetric Channel Noise (h/o = 3.16) Step edge

No filtering Averaging w-3 Median filtering w-3

Ramp edge (three points)

P,

Error size std. dev.

P,

Error size std. dev.

0.19

0.90

0.29

1.12

0.29

0.94

0.30

I .03

0.25

0.83

0.25

0.80

The analysis for the case of a ramp edge can be done in a similar manner and will not be further pursued. The same result as before is expected, however,-namely, averaging performs better than median filtering for ramp edge location estimation. We have reached the following conclusions. (i) For the two noise reduction filters we studied, the performance of leastsquares edge location estimation is degraded rather than improved if filtering is done before estimation. (ii) For step edges, average filtering is less preferable than median filtering before least-squares edge estimation. For ramp edges, the opposite is true. (Only window size 3 is explicitly discussed for median filtering). This agrees with the observation that median filter preserves sharp edges better than averaging. 9. EXPERIMENTAL

RESULTS FOR GAUSSIAN NOISE

Computer experiments support the theoretical results in Sections 2-8. The numerical results from 100 experiments are listed in Tables 2 and 3. Six typical edge sequences for location estimation are shown in Fig. 10. TABLE 5 Distribution of Error Size in Table 4 Step edge

Error size.

-4 -3 -2 -1 0 1 2 3 4 5

No fik.

Aver. filt .

Med. fik.

1 0 4 7 81 2 4 0 0 1

2 0 3 16 71 7 0 0 0 1

1 0 2 14 75 7 0 0 0 1

No Aver. fik. 3 0 4 6 71 10 5 0 0 1

Ramp &3e Med. filt.

filt.

1 2 3 8 70 11 4 0 0 1

1 0 0 10 75 23 0 0 0 1

MEDIAN

FILTERING

243

AND EDGE LOCATION

(a)

(d)

(cl

(f)

FIG. Il. Typical edge sequences in the same arrangements as Fig. 10 except that Gaussian noise is replaced by BSC noise.

244

YANG 10. BINARY

SYMMETRICAL

AND HUANG

CHANNEL

NOISEEXPERIMENTAL

RESULTS

So far, we have looked at edge sequencescontaminated by Gaussian noise. For impulse-type noise, we expect median filtering to perform better than averaging because median filtering not only preserves edges better but also reduces heavytailed (this refers to the probability density function of the noise) noise more effectively. We present the experimental results based on 100 edge sequences for each case. The noise here is due to a binary symmetric channel; each signal was transmitted by F’CM through a binary symmetric channel so that each bit has a probability Pb of being changed from 0 to 1 or vice versa. They are listed in Tables 4 and 5. One interesting observation is that median filtering improves location estimation for ramp edges but not for step edges. This is because the locations of step edges are easily shifted by 1 if one noise spike occurs around the step. In Fig. 11, we show six typical edge sequences with binary symmetrical channel noise. Because a ramp is a much better model than a step for a real-life edge, it is expected that median filtering will improve edge location estimation accuracy for real-life edges corrupted by impulse-type noise. 11. CONCLUSION

A performance evaluation of median filters was carried out by measuring their effects on one-dimensional edge location estimation (using least squares), and comparing with those of averaging. We looked at Gaussian noise (a light-tailed noise) as well as BSC noise (binary symmetrical channel noise, a heavy-tailed noise). Theoretical performance predictions were derived for the Gaussian case, and computer simulations carried out for both the Gaussian and the BSC noise cases. It was found that median filtering preserves step (sharp) edge location better than averaging, but does not preserve ramp edge location as well for the case of Gaussian noise. For BSC noise, median filtering preserves both step and ramp edge locations better. The surprising conclusion is that in the case of Gaussian noise neither averaging nor median filtering improves the accuracy of edge location estimation by least squares. In the case of impulse-type noise, median filtering does improve edge location estimation accuracy for ramp edges which are reasonable models for real-life edges. REFERENCES 1. G. J. Yang and T. S. Huang, Median Filters and Their Applicatiom to Image Processing, Tech. Rep. TR-EE 80-1, January 1980, Purdue University. 2. J. W. Tukey, Exphmatoty Data Ana&sis, preliminary ed., Addison-Wesley, Reading, Mass., 197 1 3. J. W. Tukey, Explomtoty Data AI&&, Addison-Wesley, Reading, Mass., 1977. 4. L. R. Rabiner, M. R. Sambur, and C. E. Schmidt, Applications of a nonlinear smoothing algorithm to speech prowsin g, IEEE Trans. Acowt. Spevxh Signal Process. ASSP-23, 1975,552-557. 5. N. S. Jayant, Average and median-based smoothing tecbniquca for improving digital speech quality in the presence of tiansmission errors, IEEE Tram. Comm. COM-M, 1976,1043-1045. 6. W. K. Pratt, Median filtering. in Srmionmccf Rqort, Image Processiug Institute, University of Southern California, Sept. 1975, pp. 116-123. 7. B. R. Frieden, A new restoring algorithm for the preferential enhancement of edge gradients, J. Opt. Sot. Amer. 66, 1976,280-283. 8. P. M. Namndra, “A separable median filter, Proc. Pattern Recognition and Image Proceasing, May 31-June 2, 1978.

MEDIAN

FILTERING

AND EDGE LOCATION

245

9. M. I. Shamos, Robust picture pr oossing operators and their implementation as circuits, Proc. ARPA Workshop on Image Understanding, pp. 127- 129, 1978. 10. W. L. Eversole, D. J. Mayer, F. B. Frazeo, and T. F. Cheek, Jr., Investigation of VLSI Technologies for Image Processing, Proc. ARPA Workshop on Image Understanding, pp. 191- 195, 1978. 11. T. S. Hung, G. J. Yang, and G. Y. Tang, A fast twodimensional median filtering algorithm, IEEE Trans. Acoust. Speech Signal Process. AS!W-27, 1979, 13-18. 12. S. G. Tyan, Median filters: Deterministic properties, in Two-Dimensional Signal Processing: Transfom and Median FiIters (T. S. Huang, Ed.), Springer-Verlag, Berlin/New York, 1981. 13. B. Justusson, Median filters: Statistical properties, in Two-DimemionaI Signal Processing: Transform and Median Filters (T. S. Huang, Ed.), Springer-Verlag, Berlin/New York, 1981. 14. P. F. Velleman, Robust nonlinear data smoothers: Definitions and recomm endations, Proc. Nat. Acad. Sci. US 74, 1977, 434-436. 15. J. W. Burnett and T. S. Huang, J. Opt. Sot. Amer. Image mensuration by maximum a posteriori probability estimation, Feb. 1978.