Volterra

IEEE TRANSACTIONS ON IMAGE PROCESSING, VOL. 5 , NO. 6, JUNE 1996

950

A General Framework for Quadratic Volterra Filters for Edge Enhancement Stefan Thurnhofer, Member, IEEE, and Sanjit K. Mitra, Fellow, IEEE

Abstruct- An inherent problem in most image enhancement schemes is the amplification of noise, which, due to Weber’s law, is mostly visible in the darker portions of an image. Using a special class of quadratic Volterra filters, we can adapt the enhancement process in a computationally efficient way to the local image brightness because these filters are approximately equivalent to the product of a local mean estimator and a highpass filter. We analyze and derive this subclass of quadratic Volterra filters by investigating the 1-D case first, and then we generalize the results to two dimensions. An important property of these filters is that they map sinusoidal inputs to constant outputs, which allows us to develop a new filter characterization that is more intuitive for our application than the 4-D frequency response. This description finally leads to a novel least-squares design methodology. Image enhancement results using our Volterra filters are superior to those obtained with standard linear filters, which we demonstrate both quantitatively and qualitatively.

I. INTRODUCTION HE human visual system is able to discern details in images depending on the characteristics of the surrounding area. Many different effects can influence this process, and various types of masking are described in the literature [l]. One of the most important ones is commonly known as the Weber-Fechner law, and it applies in a similar form to other sensory processes as well. It characterizes the perception of grey-level differences in dependence on the surrounding background image intensity. Roughly stated, the just noticeable difference (JND) is proportional to the average intensity of the surrounding pixels. In other words, we can see details more easily in dark regions, whereas bright portions tend to mask details. Even though Weber’s law does not describe such a dependence for very small intensity values, we can still use the proportionality as a good overall estimate. As a direct consequence of Weber’s law, image processing systems (image enhancement, image coding, image restoration, etc.) must treat the darker regions of an image very carefully because a human observer would perceive any imperfections more easily here than in the brighter areas. This is especially Manuscript received December 15, 1994; revised December 11, 1995. This work was supported by the SDIORST, managed by the Office of Naval Rescarch under contract ONR N00014-85-K-0551 and by a University o f California MICRO grant with matching funds from Hughes Aircraft Co., Xerox Corporation, and Signal Technology, Inc. S. Thurnhofer is with Lucent Technologies, Bell Laboratories, Allentown, PA 18103 USA (e-mail: stefan.thumhofer@bell-labscom). S. K. Mitra is with the Department o f Electrical and Computer Engineering, University of California, Santa Barbara, CA 93106 USA. Publisher Item Identifier S 1057-7149(96)04184-X.

I

I

edge extraction

Fig. 1. Block diagram of the unsharp masking technique.

0 -1 0

-1 -1

0 -1 0

(Laplacian I) Fig. 2.

-1 -1

-1

-1 -1

-1

-1

-1

(Laplacian 11)

1 -2 1

-2

1 -2 -2 1 (Laplacian 111)

Filter kernels for three Laplacian operators

true when images are not perfectly clean but were sent through transmission channels or were otherwise corrupted since noise in dark areas is very visible, and it should not be amplified even more. For example, consider the structure of Fig. 1. The input image is sent through a block that extracts edges and features. The output is then scaled by an appropriate factor p and added back to the original image. This method is generally referred to as unsharp masking [2],[3] and is quite effective for enhancing low contrast images. The edge extraction block in Fig. 1 is often implemented as a linear highpass filter such as a discrete Laplacian operator for which three possible filter kernels are shown in Fig. 2 [ 3 ] .The main advantage of using this particular filter is its computational simplicity. The highpass filter enhances those portions of the image that bear mostly high-frequency information, i.e., edges and textured regions. The perceptual impression is improved because the image appears sharper and better defined. An apparent problem of this technique is that it does not discriminate between actual image information and noise. Thus, noise is enhanced as well. To eliminate this problem while still preserving the simplicity of the algorithm, we make use of the property of the HVS described by Weber’s law and modify the method such that the image enhancement is dependent on the local average pixel intensity. In bright regions, we can enhance the image more because noise and other gray-level fluctuations are much less visible. On the other hand, in darker portions, we want to suppress the enhancement process since it might deteriorate image quality. This simple idea indicates the need for a highpass filter that depends on the local mean similar to H ( w ) cx H ~ ( w. )(local mean)

1057-7149/96$05.00 0 1996 IEEE

Authorized licensed use limited to: INDIAN INSTITUTE OF TECHNOLOGY KHARAGPUR. Downloaded on March 8, 2009 at 22:01 from IEEE Xplore. Restrictions apply.

(1)

95 1

THURNHOFER AND MITRA: GENERAL FRAMEWORK FOR QUADRATIC VOLTERRA FILTERS FOR EDGE ENHANCEMENT

where HtL(w)is a linear highpass filter. The most obvious solution of this problem is to multiply H h ( w ) with a linear local mean estimator (i.e., an averaging filter). The exact computational expense for such a structure depends on the particular filter design, of course. For a filter with a region of support of 3 x 3 pixels, we have to expect abjout 10 multiplications and 16 additions. We will see that using quadratic Volterra filters, a much more efficient solution can be found that only needs seven mulliplications, of which two could be implemented as shift operations and four additions. Among the various forms of nonlinear filters, Volterra filters represent the most natural extension of linear filters. Volterra filters can be described as essentially a linear filter with higher order polynomial extensions. Thus, reducing a Volterra filter to its first-order component yields a linear filter. On the other hand, this relationship with linear filters is tht. basis for some very useful properties of the higher order terms. For example, we can define the describing parameteIs as a generalized impulse response and its multidimensional Fourier transform (or Laplace transform for the continuous-time case) as a generalized frequency response. Even though the filter is not linear with respect to the input signal anymore, it is still linear in the impulse response coefficients, i.e., a linear combination of filters is equivalent to a filter with the same linear combination of the kernel parameters. Another important observation is that Volterra filters work on products of input samples. Both of these simple properties are fundamental for the analysis because it allows us to write a Volterra filter as a multidimensional convolution operation. We restrict ourselves to second-order filters only since they are the smallest-degree filter that exhibit a nonlinear behavior. Higher-order filters can potentially be very expensive in terms of number of computations, and for implementation purposes, it is desirable to keep the number of arithmetic operations as small as possible. In Section 11, we describe the basic properties and rellationships of quadratic Volterra filters. Section I11 defines various subclasses of quadratic Volterra filters, which can be vvritten as mean-weighted highpass filters in one and two dimensions. We also introduce a scheme to describe their dependence on input frequency and orientation. Based on this description, we use a least-squares approximation to design a 2-D filter for perceptual edge extraction. Finally, Section IV illustrates the application to image enhancement.

We can then write y(n) = hz(n,1,n2)* *2(n>1,n2)l n=n1=n2

where ** denotes 2-D convolution. For example, consider the system given by ~ ( n=) 2z(n - 1)x(n) - 3 ~ ( n ) x~( n l)z(n - 2 ) . Using the notation for unit-sample sequences

+

the Volterra kernel for this system is then

26(nl

36(nl,n2)

+

+ 1,1%2

A discrete-time quadratic Volterra filter is defined by 141 00

Y).(

= V2(.) =

00

hz(kl,kz)z(n-kl)z(n-k.2). kl=--oo

&=-a

+

hZ(nl,n2) =

2). Even though Volterra systems are nonlinear in the input sequence, they are linear in the kernels. It is this observation that allows us to use frequency domain techniques for the analysis and design of these systems. We define the spectrum of the Volterra kernels as the multidimensional Fourier transform of h,2(ni,n2) [SI, 161 -

1,1%2) -

6(TLl

-

nl=-w nz=-oo

The symbol .F denotes the Fourier transform in one or more dimensions, depending on the context. We call H z ( w l ,w 2 ) the second-order generalized frequency response [6] or, for simplicity, second-order frequency response. Similarly, hz(n1.n ~ is) called the second-order impulse response. Similar to the design problem of linear filters, the frequency domain representation is an important tool for the design of Volterra filters. We rewrite (4) using the frequency responses of input and system as V(n) = - T - 1 { H L ( w l . w Z ) x ( w 1 , w 2 ) } 1n=n1=n2

(8)

where X ( w ) = . F { x ( n ) }and X ( w 1 , w z ) = F { % ( n l , n z ) = } X ( w l ) X ( w 2 ) .Unless noted otherwise, we use lowercase letters for sequences in time or space and uppercase letters for their Fourier transforms. Unlike for linear filters, however, the impulse response of quadratic Volterra filters does not necessarily define the filter uniquely, i.e., there can exist many kernels for the same filter [7]. Only if the kernel is symmetric does it uniquely describe the system, and for a second-order system, this means that [8] G(n1, n2) = q n 2 , n l )

11. QUADRATIC VOLTERRA FILTERS: DESCRIPTION AND RELATIONSHIPS

(4)

(9)

which we will use frequently. Since every kernel can be symmetrized [SI, (9) does not impose any restrictions. The 2-ID Volterra series is a straightforward extension of (2) [9]. Here, the input and output signals are 2-D sequences (usually called “images”). We write the 2-D discrete-time quadratic Volterra filter as

(2) This expression can be interpreted a.s a 2-D convolution of hz(n1,n2) with products of samples of the input sequence. This becomes clearer if we define

Authorized licensed use limited to: INDIAN INSTITUTE OF TECHNOLOGY KHARAGPUR. Downloaded on March 8, 2009 at 22:01 from IEEE Xplore. Restrictions apply.

IEEE TRANSACTIONS ON IMAGE PROCESSING, VOL. 5, NO. 6, JUNE 1996

952

with 3(nJ1,n2,n3,nq) = z(n1,n ~ ) z ( n n 3 g, ) . We assume that symbols like zor y can represent both 1-D and 2-D sequences, depending on the context. In (lo), we have omitted the limits of the sums. In the most general case, all variables kl,. . . ICg will range from -cc to 00. All the filters that we consider here are of finite extent, however. If, for instance, the impulse response of the second-order term is such that the filter only uses pixels that are a maximum of N I pixels away from the current pixel in n,l direction and N2 pixels away in n2 direction, then - N I 5 k l , k3 I N I , and -Nz 5 k z , kq 5 N2. In general, the limits are of minor importance, and in the following, we will omit them for simplicity unless their values cannot be inferred easily. Using again the Fourier transform notation, we can rewrite (11) as

if

The approximation error is then given by

where X ( w ) = X ( W ) - 27rpL,S(w) is the mean-removed version of X ( w ) . Pro08 The proof of this theorem is given in [12]. The basic idea is to decompose the general expression of quadratic Volterra filters given in (2) into the two terms

111. MEAN-WEIGHTED HIGHPASSFILTERS

We now define a subclass of quadratic Volterra filters and derive its properties. The filters from this class are approximately equal to a product of a local mean estimator and a highpass filter. We investigate the 1-D case first since the arithmetic expressions are simpler than those for the 2-D case. A. I - D Mean- Weighted Highpass Filters

In [lo], a special quadratic Volterra filter called Teager's algorithm has been defined, and in [11], it has been shown that this filter is approximately equal to a product of local mean and highpass y(n) = ."n)

-

k,(Zz(n)

+

l)z(n z ( n - 1)

z(n -

-

+

+ 1) -

+

z ( n + 1))

(13) (14)

where k , = ( z ( n- I) z ( n ) x(n 1))/3 estimates the local mean. Additionally, 2-D extensions of the Teager filter have been introduced in [ l l ] , which were found intuitively and could be decomposed analogous to (14). These filters have been used successfully for image enhancement. However, in [l I], it was not clear whether or not these filters would be optimal in any way or what is the general class of Volterra filters to which they belong. We can now answer the latter question for 1-D filters with the following theorem, where we model the input sequence z ( n ) by an i.i.d. process with uniform distribution and with mean px and standard deviation a,. Theorem: A 1-D second-order Volterra filter can be approximated by a mean-weighted filter, i.e.

ki

ka

where ? ( n )= z ( n )- p x is the mean-removed version of z(Tr1). We then show under which conditions yz(n) can be neglected 0 compared with yl(n). We need to make some remarks about this theorem. First, (16) is true for practically all images. In those cases where it should not be true, the signal or image would have to be scaled and shifted accordingly before processing. Second, we can approximate the mean p, in (15) by a local averaging filter if we assume that the kernel of the Volterra filter is finite. In this case, only pixels in a local neighborhood influence the result, and thus, only the local characteristics of the signal play a role for the computation of the output. Third, there is no restriction on the spectral shape of Hz(w,0). For applications in edge enhancement, however, we will require highpass characteristics, in which case, the corresponding Volterra filter is a mean-weighted highpass filter. For the most general case, we can relax the conditions in (16) and (19) even further and drop the symmetry condition in (18), but then, we need to rewrite (15) as

Even though (19) appears fairly complicated, we can greatly simplify it for certain special cases. In the first class, filters consist of squares of pixels. Dejnition: If h~(n1,na) = 0 for all n1 # nz and E,, hz(n1,nz) = 0 and Hz(w,O)has highpass

E,,

Authorized licensed use limited to: INDIAN INSTITUTE OF TECHNOLOGY KHARAGPUR. Downloaded on March 8, 2009 at 22:01 from IEEE Xplore. Restrictions apply.

953

THURNHOFER AND MITRA: GENERAL FRAMEWORK FOR QUADRP,TIC VOLTERRA FILTERS FOR EDGE ENHANCEMENT

and thus

The locations of the nonzero coefficients of the impulse response are indicated in Fig. 3(b). These filters have the additional useful property that they are symmetric about the current sample. The symmetrized Teager filter y(n) = x 2 ( n )0.5z(n- l ) z ( n 1)- 0.5x(n 1)x(n- 1) is an example for a simple system from this class [12]. The Teager filter maps sinusoidal inputs to constant outputs [lo], and it is straightforward to show that, in fact, every class I1 system has this property. Using (21), we can see that H z ( w , w ) = C h 2 ( k , - k ) = H 2 ( 0 , 0 ) , which is assumed to be zero. With X ( w ) = -j7rS(w - W O ) j n S ( w wg) and the spectrum of the output signal given by [12]

+

(a)

(b)

Fig. 3. Locations of nonzero impulse response coefficients of special cases of mean-weighted highpass filters: (a) Class I filters; (b) class I1 filter.

+

+

characteristics, then the filter is called a class I discrete-time quadratic Volterra filter, or class I filter for short. The following lemma states the approximate behalvior of class I filters. Lemma I : Class I filters can be approximated as ZL meanweighted highpass filter. Proof: Equation (17) holds by assumption, and (18) is true since we have hz(n1,na) = Cy=-, h 2 ( k , k ) 6 ( n l k , n 2 - k ) , which leads to 00

Hz(w1, w2) =

h 2 ( k ,k ) e - + - J k W z

+

the output consists only of a dc component.

B. 2 - 0 Mean- Weighted Highpass Filters We now consider the extension of the previous results to 2D signals for applications in image processing. Since the class I1 filters of Lemma 2 inherently lead to symmetric systems, the localization of edges is preserved when we apply class I1

f=-00

filters to images. This property is especially important here because otherwise, it would be impossible to enhance edges and thus H z ( w , 0) = h z ( k , k ) e - J k w = Hz(0, w ) . It remains to show that S h 2 0. In (19), all the sums since the exact location of the feature would be unknown. The fundamental building block of class I1 systems is yield zero because h z ( k l , k2) is nonzero only for kl = k2. z(n - k)x(n k ) . The overall system consists of a linear Therefore, we obtain Sh = 0 here. 0 Examples for this class are y(n) = z 2 ( n )- 0.5z2(n - 1) - combination of several of these terms for different k . Thus, 0.5x2(n-3) and y(n) = z 2 ( n - 1 ) - - x 2 ( n - 2 ) . The first filter each term is a product of samples that are equally far away to can approximately be written as p z ( 2 x ( n )- x ( n - 1)-- z ( n- the left and right of the center location n. We extend this idea 3)) and the second as 2 p z ( x ( n - 1) - x(n - 2)). Fig. 3(a) to the 2-D case simply by defining the basic building block of shows the possible locations for nonzero coefficients of the the 2-D filter as z(n1- k l , nz - k z ) z ( n l + k l , 722 k z ) , i.e., it consists of pixels centered around the current pixel z(n1,n2). impulse response in the nl-n2-plane. The second family of filters is more important for our For the impulse response, we obtain the simple property applications. h , 2 ( n l l n 2 , n 3 , n 4 )# 0 only if n1 = -n3 and n2 = -724. Definition: If h 2 ( n l r n 2 ) = 0 for all n1 # -712 and (23) E,, E,, h2(a1, 7L2) = 0 and H2(w, 0) has highpass char- Analogous to Lemma 2, we also require that acteristics, then the filter is called a class I1 discrete-time quadratic Volterra filtei, or class I1 filter for short. The following lemma states the approximate behalvior of k.1 k.2 k3 f4 class I1 filters. Additionally, we can show that Lemma 2: A Class I1 filter can be approximated as a meanweighted highpass filter. Proof: Here, we can again simplify S h , and with the To prove this, we write the kemel as symmetry condition for hz ( k l , k z ) , we obtain

+

+

S" = x h ; ( k , - k ) 2 0. k k#O

In a similar fashion as in the proof of Lemma 1, we can show that

Then, we find that

00

k=l

Authorized licensed use limited to: INDIAN INSTITUTE OF TECHNOLOGY KHARAGPUR. Downloaded on March 8, 2009 at 22:01 from IEEE Xplore. Restrictions apply.

954

IEEE TRANSACTIONS ON IMAGE PROCESSING, VOL. 5 , NO. 6, JUNE 1996

and thus

a=- 1 a
a
where we have used the symmetry of the kernel, i.e.

t

ko Using (23)-(25), we find the corresponding extension to (15) to be

has highpass characteristics in Assuming that Hz(w1. w2,0,0) either the w1 or w2 direction or in both, we conclude that the class I1 2-D systems can be approximated by mean-weighted 2-D highpass filters. C. Characterization of Edge Extracting Filters

We have developed the theory behind the class of filters that we want to use for edge extraction. We know their structure and properties in both space and frequency domains, and we also have justified why, at least in principle, this family would be advantageous for our goal. In the next step, we need to develop a performance measure that somehow allows us 1.0 find the best filter from this class. Even though essentially all mean weighted highpass filters would extract edges of an image in some way, there are still important differences that we have to take into account. For instance, a problem that would not occur in I-D signal processing but is of extreme importance for image processing is isotropy. Of course, the filter should find edges independent of their orientation. Horizontal or viertical boundaries must lead to the same response, and therefore, the degree to which a filter is isotropic is critical for the resulting image quality. In this section, we introduce a method of characterizing the isotropical behavior as well as the dependence of a filter output on the input frequency [13]. Unlike for linear filters, however, input frequencies have a rather complicated influence on the output, and we cannot show this dependence in any cornplete and yet easily interpretable way. Even though the filter is mathematically completely defined by the frequency response H2 (.), a plot of this function does not usually yield any insight into the edge extracting capabilities of the system. Essen tially, two factors make the analysis hard to interpret: the complexity of the system and the complexity of the input signal. Even for simple quadratic systems, the frequency domain representation can appear confusing. Above all, we cannot restrict our design to systems for which we have a simple intuitive explanation. This leaves us with the second source of complexity: the input signal. Of course, we have complete control over this for design purposes, which means that we can develop a filter that reacts in a specific way to a certain input. However, we have to make sure that the filter's response to real-world signals

a=l b=O

b>O

Fig. 4. Parameters a and b determine the orientation of the input sinusoid. They are interlinked by (I' h' = 1.

+

(i.e., arbitrary images) remains predictable. This simplification allows us to describe the behavior in more intuitive terms than by using the complete frequency response. The choice of the input is very much determined by our design goals. The signal must be such that the system response is meaningful so that we can draw conclusions as to the extent to which we have achieved the desired behavior. To this end, we choose a single sinusoidal signal as the system input. By varying its frequency, we can roughly estimate how the system will respond to more complex input signals with arbitrary frequency content. Even though it will not be possible to predict quantitatively how a system will react to an arbitrary input signal, because we do not take into consideration any intermodulation frequencies or harmonics, we still can make qualitative assertions. If the region of support of the filter is small enough, we can model the corresponding portion of the input image by a single 2-D sinusoid. The approximation error will decrease with the size of the region of support (typically 5 x 5 or 3 x 3). As we have already indicated at the beginning of this section, isotropical behavior is of critical importance. To incorporate a test for isotropy, we use rotated sinusoids with the rotation as a free parameter

+

where a and b determine the orientation, and ' a b2 = 1. In Fig. 4, the possible orientations are indicated with the corresponding values for these parameters. Both must be between -1 and 1, and because they are interlinked, only one of them is sufficient to define the orientation. In addition, for our purposes, half of the possible directions in Fig. 4 are redundant, and we will only use

h = 2/1 - a 2 . The output spectrum of a 2-D quadratic Volterra filter is given by [12]

Authorized licensed use limited to: INDIAN INSTITUTE OF TECHNOLOGY KHARAGPUR. Downloaded on March 8, 2009 at 22:01 from IEEE Xplore. Restrictions apply.

955

THURNHOFER AND MITRA: GENERAL FRAMEWORK FOR QUADRATIC VOLTERRA FILTERS FOR EDGE ENHANCEMENT

The corresponding frequency responses are

We find

H;a)(W1,W2,Wj,w4) = (1 - cos(w1

- Wj))

x (1 - cos(w2

- wg))

(35)

Y y ( a ) ( w O ,U ) = 2 7 4 1 - cos(woa))

x (1 - cos(wob)).

(36)

Fig. 5(a) shows the isotropy plot. The right-hand plot displays the same graph viewed from the side. Here, the U axis is no longer visible since it is perpendicular to the paper plane with the positive part pointing into the paper. This side view makes As expected, the output consists of components with zero and the dispersion clearly visible, i.e., we can assess the extent to twice the input frequency. We can, however, easily show that which the filter varies its response for different values of U . only the dc component will be nonzero. We use (26) and (23) We see that Filter A does not respond at all for 0, = -1, and obtain 0 and 1, corresponding to horizontal and vertical orientations (see Fig. 4). Similarly, the image in Fig. 6(a) shows the strong -h, - k z ) == 0. directional dependence of the filter output. We can make another observation by examining the way the response rises for increasing frequencies in Fig. 5(a). For small values of Thus, only the terms with S(wl,wz) remain in (31). This is W O , the output is very small, and it starts to rise only for W O not too surprising because the class I1 2-D filters were derived approximately greater than 0.5 rad. This means that the system as an extension of their 1-D counterparts in Section 111-A and has a strong highpass characteristic. Accordingly, the lines in for those we have already shown this property. Thus, we write Fig. 6(a) are very thin. In [ 111, the authors introduced heuristically found extensions of Teager's algorithm. Now, we can study their properties in the framework of the isotropy plots. The first with ~i'~)(wg, a ) = 7 r 2 [ ~ 2 ( w o aweb, , - q a , -web) of these systems (Filter B) is an application of the Teager Hz(-woa, -web, w o ~ , , w o b ) ] . filter along diagonals. Summarizing, we have developed a method that allows us db)(71.1, 71.2) to characterize the behavior of class I1 filters in an intuitive = 2zy71.1, 712) - z(n1 1.712 1) way. Only three values play a role in this scheme: the input frequency, the orientation of the input sinusoid, and the level x z(n1 - 1 , 7 1 2 - 1) - z(n1 1,712 - 1) of the constant output. Most importantly, we can show the x x(n1 - 1,712 1) (37) relationships graphically, which is a fundamental step toward H;"@l, w2. wx,w4) intuitive understanding. Our assumptions greatly simplify the originally rather complex situation. = 2 - cos(w1+ w2 - wj - w4) We now want to illustrate how to interpret the plots of - cos(w1 - wz - wg wg) (38) directional and frequency dependence. We call these graphs Y;c(b) (WO, a) isotropy plots because they indicate to which degree i2 filter is isotropical. In order to understand how these properties = 2 7 4 2 - cos(2woa 2 w o J I - a 2 ) determine the system output for images, we also use an - COS(2W"U - 2w02/1 - U * ) ) . (39) artificial image shown in Fig. 6(a). It is a ring with added Fig. 5(b) shows that this system is fairly isotropic up to about Gaussian noise ( 0 2 = 0.01). The image dimensions are 256 x 256 pixels, and the inner and outer radii of the ring are 80 wo = 0.5, which we can verify visually in Fig. 6(b). Filter C is similar to Filter B, but along the horizontal and and 100 pixels, respectively. Before adding the noise, pixels on vertical directions, and it is defined by the ring have unity value, and outside, they are zero. The step edge is almost ideal, thus containing very high frequencies, y(C)(n1, n 2 ) 2x2(711,n z ) - z(n,1 1,n 2 ) which will illustrate the differences between the filters more x z(71.1 - 1,'TLz) - Z(71.117L2- 1) clearly than a smooth transition. x 5(7L1. r12 1) (40) We define the first example system (Filter A) by the inputoutput relation H$')(wl,w2,"-'3, '"4)

+

+

+ +

+

+

+

+

+

= 2 - COS(W~- ~

Y;c(c)(wo, U)

g -) C O S ( W Z - wg)

(41)

= 27r2(2 - cos(2woa)

cos(2wo2/1 - a").

(42)

The isotropy plot in Fig. 5(c) appears even smoother than for Filter B (see (37)) and shows isotropical behavior up to higher frequencies. In addition, Filter C has a more pronounced

Authorized licensed use limited to: INDIAN INSTITUTE OF TECHNOLOGY KHARAGPUR. Downloaded on March 8, 2009 at 22:01 from IEEE Xplore. Restrictions apply.

956

IEEE TRANSACTIONS ON IMAGE PROCESSING, VOL. 5, NO. 6, JUNE 1996

w .......

3

..; .... ....................... . . . . . .

~

25

a

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

:.

. . .

............. ..-.

a

a

Fig. 5. Isotropy plots for the example filters. The left graph shows the regular view, and the right one displays it from the side with the a-axis perpendicular the paper plane.

highpass characteristics, which explains why the lines in Fig. 6(b) are thicker than in Fig. 6(c). Summarizing, we find that for the design of an edge extracting filter, the first system would not be a good choice because of its very unisotropical properties. Either Filter B or C would be better options, where Filter C is slightly more sensitive to details.

D. Least-Squares Design

ckl

with Ck2hz(k.1,k 2 , - k l , -kg) = 0. Obviously, any filter from this class is a linear combination of terms like z(nl - k l , ri2 - ~ Z ) Z ( T L ~ kl,ng k2). This is based on the property of Volterra filters that they are linear in the kemels. We make use of this observation and introduce the concept of basis filters for the design of Volterra filters. We define an expression

+

yicl

k 2 (711.712)

By combining (10) with (23), we obtain the input-output relation for class I1 filters [14] Y(”1.122)

= kz

z ( n l - Icl) n2

- k2)x(nl

-

x(n1

x z(n1 - k l , rL2

-

+ kl, + k z ) n2

k2)

+ ICl n2 + k 2 ) 1

(43)

(44)

as a basis for class I1 filters, where 0 5 k1 < m, -00 < k 2 < 00, and k l 1k.21 # 0, i.e., kl and kg cannot be equal to zero at the same time. We only consider unique pairs of kl and kz. If we would also permit ICl < 0, then some of the filters would be identical because ykl,lc2(n1,n g )= y--lcl,-kz(n1,n2). Using

+

J L Z ( h , k . 2 , -k1, - k . 2 ) kl

= x2(n1,112)

+

Authorized licensed use limited to: INDIAN INSTITUTE OF TECHNOLOGY KHARAGPUR. Downloaded on March 8, 2009 at 22:01 from IEEE Xplore. Restrictions apply.

957

THURNHOFER AND MITRA: GENERAL FRAMEWORK FOR QUADRATIC VOLTERRA FILTERS FOR EDGE ENHANCEMENT

Fig. 6.

Original “ring” image and results of applying the example filters to it.

(44), we find the equivalent expression to (43) as Y ( n l > n 2 )= ~ ~ B k l , k Z / B l , k 2 ( ~ l , ~ 2 (45) ) kl

kz

where the O k l , k z represent the design coefficients. We prove this equivalence by rewriting (43) as

this set. Furthermore, the basis filters have the advantage that for each of them the sum of their coefficients equals zero, and thus, (17) holds by design, and the coefficients i9kl , k 2 can have arbitrary values. The relation between the kernel coefficients and the t?kl,kL is

7;

h z ( O , O , 0,O) = JL2(n1,nz1-n1,

(46)

k , k 2

ki

kz

-4 = - @ n l , n z for

nl

+

(47) In21

#0

O
ki

kz

ki

k2

with the restrictions for 5 1 and k 2 as before for (44)1, and where we have substituted Okl,kz = - h z ( k l , k 2 , - k l , - 5 2 ) . Therefore, the set of filters in (44) completely describes all possible filters from class 11, and it is not redundant since none of the filters in (44) can be written in terms of any other from

-n2)

=0

for

nl

-oo
< 0.

(48)

In the frequency domain, the linear dependence on the basis filters is preserved by the Fourier transform. Since we know the frequency characteristics of each of the basis filters, we must find the coefficients in such a way that the overall system has a certain desired frequency response. We need to study the basis filters in more detail. We are interested in designing class I1 systems that operate on a relatively small neighborhood. This keeps the computationally complexity low for both design and implementation. We choose a 5 x 5 region of support because it is a reasonable compromise between the computational complexity of the filter and the degree of freedom for the design. We also have to make sure that modeling the input signal by a single sinusoid

Authorized licensed use limited to: INDIAN INSTITUTE OF TECHNOLOGY KHARAGPUR. Downloaded on March 8, 2009 at 22:01 from IEEE Xplore. Restrictions apply.

9%

IEEE TRANSACTIONS ON IMAGE PROCESSING, VOL. 5 , NO. 6, JUNE 1996

remains a valid characterization and, thus, should keep the region of support as small as possible. We obtain a total of 12 basis filters for -2 5 kp 5 2 and 0 5 k1 5 2 in (44). We find ~ 1 ,ws, w g ) as their frequency responses H i k 1 ' k 2 ) ( w2,

H i 0 q w l ,wp>wI(,w g )

=

H ; o ' 2 ) ( W 1 , w2, W Q ,w g )

= 1 - cos(2w2

1 - COS(W*

-

wg)

-

24)

H$l,-q)(wl, w.2;w 3 , w g ) = 1 - cos(w1

-

2w2

H;l 1 - 1) ( W l , wp, wg; w g )

-

w2 - wg

1 - COS(W1

=

H p o ) ( w l ,w 2 , W Q ,w q ) = 1

-

H ; 1 , 2 ) ( w lw.2; , w3, wg) = 1 1-

2)

( w 1 ,w.2;wg, w g )

H . p ( w 1 w.2

~

wg)

w3, w g )

-

-

-

cos(2w1

= 1-

+2 4 )

+

wg)

-

wg)

-

w3 - 2 w 4 )

-

-

2

+ 2wg) ~ + 3 wg)

2w3

2w3)

+ cos( 2w1 + 2wz

= 1 - cos(2w1

w3

-

2w2

1 - C O S ( ~ W-~ wp

H;2,")(w1:wp;w 3 , w q ) = 1 H ' p ) ( W I wp; WJ,

+ wp cos(w1 + 2w2

= 1 - COS(2Wl

- 1)( ~ 1~,p , w y , w g= )

w3

- w3)

1 - cos(w1

H p ) ( w l , w z ; W j , W q )=

H;2

cos(w1

-

wq -

2w3

-

-

2w3

wq)

-

TABLE I

k ~

0.1 0 0.01323 1.01641 1 0.00061 0.04733 2 -0.00011 -0.00897 0.00616 0.47368 3 0.01325 1.01759 4 5 0.00616 0.47368 6 -0.0001 1 -0.00897 7 -0.00021 -0.0 1674 8 -0.00012 -0.00961 9 0.00054 0.04197 10 -0.00012 -0.00961 11 -0.00012 -0.01 674 OT

0". 1.o 0.0 0.0 0.5 1.o 0.5 0.0 0.0 0.0 0.0 0.0 0.0

We write the linear combination of the basis filters with coefficients Q k l , k a as a matrix product BO, where

2 4 .

We use the isotropy plots that we developed in Section 111-C as a characterization of class I1 filters. For each of

the basis filters, we obtain the response to a sinusoidal input. For simplicity, we denote them B k l r k E( W O . U ) instead of y;C,(kl,ka)

(WO

, a).

Thus, we can express the design problem as the linear matrix equation

z=BH+e.

Bo,1(wo, U) = 27r2(1 - cos(2wob)) B",Z(WO, U ) =

27r(1

-

cos(4wob))

B1,-2(wo, U ) =

2 7 4 1 - cos(2woa

-

4wob))

131,-1(Wo, U ) =

27r(l

-

2wob))

-

cos(2woa

B1,o(wo, U) = 2 7 r ( 1

-

cos(2w"a))

2n2(1

-

COS(2W"U

B1,1(wo, U) =

B1,2(wo, U) = 2 x 7 1 - cos(2woa

+ 2wob))

+ 4wob))

B~,-~(wO, U)

= 27r2(1 - cos(4woa - 4wob))

B2,-1(wo, U )

= 2 7 4 1 - cos(4woa

(51)

The desired behavior of the ideal filter is contained in z, which we sample in the same fashion as described above for B k l k z ( w ( ) .a ) , and e represents the error vector. We find the optimal solution 8, which minimizes the error in the leastsquares sense [lS]

H^ = [B'B]

-'BTz.

(52)

We want to design a 2-D version of the Teager filter. In [lo], it has been shown that Teager's algorithm responds with an output y ( n ) = sin2(wo) if the input is a sinusoid B~,o(w U) o= , 27rz(1 cos(4woa)) z ( n ) = sin(won). Now, we need to find a 2-D system that B ~ , l ( w oa,) = 27r2(1 cos(4woa 2 w o b ) ) outputs a constant signal y(n1, na) = sin2(wo) for a sinusoidal B ~ , ~ (Uw) = o ,2?(1 - c o s ( 4 w o ~ 4wob)) input and is, of course, independent of the orientation of the excitation. where b = d m .Then, we sample these continuous We set z ( w o , u ) = sin2(wo) and convert z ( w o , a ) into the functions in -1 5 U 5 1 and 0 5 W O 5 %. Since images column vector z. Using (52), yields the parameter vector QT tend to bc oversampled, very little information is contained given in Table I. If we relax the restrictions somewhat and in the very high frequencies, and for this design, we only allow the filter output to be scaled, i.e., z(wg; U) = aisin2(wo), consider frequencies up to We choose a sample spacing of we can normalize the parameters such that E, 8 ( i ) = 3. This Aa = 0.05 and Awo = which yields 201 . 26 samples. yields 8, in Table I. Fig. 7(a) shows the isotropy plot for We order the samples into column vectors in the following manner. The first 201 elements belong to W O = 0 and from this filter. The design goal of virtually perfectly isotropical behavior has been achieved up to W O = 5. If we had used the U = -1 up to a = 1. After that, we take the 201 samples for entire range from 0 to 7r for the optimization, then we would W O = & and so on until we reach W O = and a = 1. We have obtained a less accurate behavior for smaller frequencies. number the functions B k l . k 2 ( W O , U ) in the order given before The advantage of normalizing the parameter vector is that from 1 to 12 and obtain vectors bl through b 1 2 , which we the elements of 8, can easily be approximated by simple combine into the matrix B numbers. The right-most column in Table I shows the approxB = (bl b p . .' b i z ) . (49) imation 4, using only the coefficients 1 and 0.5. Thus, we -

2wob))

-

-

+ +

5.

5,

Authorized licensed use limited to: INDIAN INSTITUTE OF TECHNOLOGY KHARAGPUR. Downloaded on March 8, 2009 at 22:01 from IEEE Xplore. Restrictions apply.

THURNHOFER AND MITRA: GENERAL FRAMEWORK FOR QUADRATIC VOLTERRA FILTERS FOR EDGE ENHANCEMENT

959

a

a

Fig. 7. Ikotropy plots for the (a) 2-D Teager filter and (b) its approximation.

obtain the input-output relationship for the approximate 2-D Teager filter

This equation is considerably simpler to implement than the full description in 6,. At the same time, the behavior of this system as shown in Fig. 7(b) does not deviate too much from the ideal one. It is still extremely isotropic for almost the entire range of interest from 0 to 7r/2. I v . IMAGE PROCESSING APPLICATIONS In this section, we want to demonstrate the application of quadratic filters to the processing of real images. The designs in the previous sections were based on the idealized assumption that the input signal is a single sinusoid, which, of course, is not true for a real image. The small region of support of the filter, however, makes this approximation reasonable if we use the instantaneous frequency for the input sinusoid. The test image is a portion of “Lena” in Fig. 8(a) and (b) shows the output of the filter in (53). Details in the brighter parts of the image have been extracted, but in the darker areas,

the filter output is very small. Accordingly, undesired noise enhancement only takes place in the brighter portions, where it is not as perceivable due to Weber’s law. The result of enhancing the image with the unsharp masking algorithm in Fig. I, and the quadratic Volterra filter in (53) as the edge extraction filter is shown in Fig. 8(c). This image is noticeably sharper, and even small details in the hair and feather regions have been enhanced. Another example is shown in Fig. 9. Again, the result of unsharp masking with our Volterra filter in Fig. 9(c) is noticeably sharper. Details have been extracted proportionally to the local image brightness. In the darker portions of the images (e.g., the tree region), the filter output in Fig. 9(bj is relatively small, whereas in the brighter portions, even minute details result in a strong response. In order to compare these results quantitatively with a linear enhancement scheme using one of the commonly used Laplacian filters, we need a criterion for image sharpness. We have found that the mean magnitude of the gradient image (using Sobel’s operator) corresponds well with the perceived sharpness. In other words, if two images appear to be equally sharp, their mean gradient magnitude is approximately the same. The measure that we use to assess the visibility of noise is a relative one, that is, it needs a reference image for comparison. It observes the portion of the local grey-level fluctuation that is already present in the reference (original) image and the part that is due to noise from the enhancement process. It also considers noise masking by the surrounding area of the image.

Authorized licensed use limited to: INDIAN INSTITUTE OF TECHNOLOGY KHARAGPUR. Downloaded on March 8, 2009 at 22:01 from IEEE Xplore. Restrictions apply.

IEEE TRANSACTIONS ON IMAGE PROCESSING, VOL. 5 , NO. 6, JUNE 1996

960

(c)

Fig. 8. Original image “Lena” (a) and edge image obtained by applying the quadratic Volterra filter in (53) to the image in (a). Note that in darker image areas, fewer edges are extracted than in the bright ones, which can be explained by the mean-weighted highpass characteristics of the filter. The result of the unsharp-masking enhancement is shown in (c) with p = 0.002.

In highly active areas, the human visual system perceives noise to a lesser degree than in flat ones. The local average luminance influences the visibility as well. Altogether, we combine these effects into a single number

(54) The index 1 stands for local and indicates that the variances and the mean are dependent on n1 and 712 and must be computed over a local neighborhood (5 x 5 pixels for our simulations). The sum is carried out over all pixels in the image. Even though the human eye cannot adapt to very small portions of the image, we assume that the observer would view the image at arbitrary distances, which justifies the use of a localized measure. The reference (original) image is denoted

by x and the enhanced one by y. Note that a comparison of the original image with itself yields k~ = 0. The value of k~ increases with the degree of image degradation and is (at least theoretically) not bounded. Since the original “Lena” and “Park” images are very clean, white Gaussian noise ( p n = 0.0, 02 = 100.0) was added to both of them. We then enhanced these images using the Volterra filter in (53), the Volterra filter C in (40) (which is called Type 1A in [ l l ] ) , and the Laplacian filters, such that the resulting overall image sharpness is the same in all corresponding cases. Table I1 lists the quantitative comparisons. Adding noise to an image does not increase pllg1l significantly, which corresponds well with our visual observations. We then computed the noise visibility measure k,v relative to the noisy image before enhancement. The

Authorized licensed use limited to: INDIAN INSTITUTE OF TECHNOLOGY KHARAGPUR. Downloaded on March 8, 2009 at 22:01 from IEEE Xplore. Restrictions apply.

~

THURNHOFER AND MITRA GENERAL FRAMEWORK FOR QUADRATIC VOLTERRA FILTERS FOR EDGE ENHANCEMENT

96 1

(c) Fig. 9. (a) Original image “Park” and (b) edge image obtained by applying the quadratic Volterra filter in (53) to the image in (a). The result of the unsharp-masking enhancement is shown in (c) with p = 0.002.

Volterra filter C shows only slightly increased noisyness since it also belongs to class of mean-weighted highpass filters. The Laplacian filters, however, lead to values of k~ that are considerably larger than for the Volterra filter in (53), which corresponds to a significant increase in noisyness. To show an example of this problem, Fig. 10 compares the results for the “Park” image visually. The result using the Volterra filter C is similar to the image in (d). Even though the three enhanced images have the same overall sharpness, we can see that Figs. 10(b) and (c) are more noisy than (d), which is most noticeable in the dark portions of the trees.

TABLE 11

“Lena”

“Park”

V. CONCLUSION We have developed a subclass of 1-D quadratic Volterra filters that can be written as a product of local mean and linear highpass. Even though some filters of this class have been

found before heuristically, we are now able to develop them in a unified framework. We have also extended this concept

Authorized licensed use limited to: INDIAN INSTITUTE OF TECHNOLOGY KHARAGPUR. Downloaded on March 8, 2009 at 22:01 from IEEE Xplore. Restrictions apply.

IEEE TRANSACTIONS ON IMAGE PROCESSING, VOL. 5, NO. 6, JUNE 1996

962

(d)

(c)

Fig. 10. Comparison of noise enhancement using u n s h q masking. The “Park” image with added Gaussian noise is shown in (a). The results of unsharp masking using the first Laplacian filter in (b) ( p = 0.45),the second Laplacian filter in (c) ( p = 0.17),and the Volterra filter in (53) in (d) ( p = 0.002) demonstrate that linear filters amplify noise, especially in the dark areas more than our Volterra filters. All the enhanced images have approximately the same sharpness. The Laplacian I11 filter performs very poorly, and the resulting image in not shown here.

to 2-D systems obtaining the general class of filters for edge extraction. Describing these systems using the kemel and the frequency response alone, however, does not enable us to find the coefficients in a simple and robust way. Thus, we have developed an approximate description that is more suitable for the problem at hand by using simple 2-D inputs. This also allows us to evaluate the orientational dependence of the filter. We have designed an optimal filter using basis filters and a least-squares frequency domain approach. Compared with other methods [4], [ 6 ] , [9], [16], [17], this technique gives a more intuitive understanding of the filter behavior and does not depend on heuristic assumption for the impulse response. As an immediate application in image processing, we have considered the unsharp masking enhancement scheme. Using

quantitative criteria, we have shown that quadratic Volterra filters yield perceptually better results than standard linear filters.

REFERENCES [ 11 K. R. Boff, L. Kaufman, and J. P. Thomas, Handbook of Perception and

[2]

131 [4]

[5]

Human Performance, vol. I: Sensory Processes and Perception. New York: Wiley, 1986. J. S. Lim, 2 - 0 Signal and Image Processing. Englewood Cliffs, NJ: Prentice-Hall, 1990. A. K. Jain, Fundamentals of Digital Image Processing, Englewood Cliffs, NJ: Prentice-Hall, 1989. G. L. Sicuranza, “Quadratic filters for signal processing,” Proc. IEEE, vol. 80, pp. 1263-1285, Aug. 1992. I. Pitas and A. N. Venetsanopoulos, Nonlinear Digital Filters. Boston, MA: Kluwer, 1990.

Authorized licensed use limited to: INDIAN INSTITUTE OF TECHNOLOGY KHARAGPUR. Downloaded on March 8, 2009 at 22:01 from IEEE Xplore. Restrictions apply.

THURNHOFER AND MITRA: GENERAL FRAMEWORK FOR QUADRATIC VOLTERRA FILTERS FOR EDGE ENHANCEMENT

161 S. A. Billings and K. M. Tsang, “Spectral analysis for nonlinear systems, Part I: Parametric nonlinear spectral analysis,” Mech. Syst. Signal Processing, vol. 3, no. 4, pp. 319-339, 1989. [7] L. 0. Chua and C. Y. Ng, “Frequency domain analysis of nonlinear systems: general theory,” Electron. Circuits Syst., vol. 3, pp. 165-185, July 1979. [XI M. Schetzen, The Volterra and Wiener Theories of Nonlinear Systems. New York: Wiley, 1980. [9] G. F. Ramponi, G. L. Sicuranza, and W. Ukovich, “A computational method for the design of 2-D Volterra filters,” IEEE Trans. Circuits Syst., vol. 35, pp. 1095-1102, Sept. 1988. [IO] J. F. Kaiser, “On a simple algorithm to calculate the ‘energy’ of a signal,” in Proc. IEEE Int. Con$ Acoust., Speech Signal Processing, Albuquerque, NM, Apr. 1990, pp. 381-384. [ l I] S. K. Mitra, H. Li, I. S. Lin, and T.-H. Yu, “A new class of nonlinear filters for image enhancement,” in Proc. IEEE Int. Conf: Acoust, Speech Signal Processing, Toronto, Ont., Canada, 1991, pp. 252552528, 1121 S . Thumhofer, “Quadratic Volterra filters for edge enhancement and their applications in image processing,” Ph.D. dissertation, Univ. of Califomia, Santa Barbara, 1995. [I31 S. Thumhofer and S. K. Mitra, “Volterra filters for perceptual edge extraction,” in Proc. 28th Ann. Asilomar Confi Signals, Syst. Comput., Pacific Grove, CA, Nov. 1994, pp. 736-740. [14] S. Thurnhofer and S. K. Mitra, “Designing quadratic Volterra filters for nonlinear edge enhancement,” in Proc. Inl. Con$ Digital Signal Processing, Limassol, Cyprus, June 1995, pp. 320-325. [ 151 J. M. Mendel, Lessons in Digital Estimation Theory. Englewood Cliffs, NJ: Prentice-Hall, 1987. [ I 61 G. Ramponi, “Bi-impulse response design of isotropic quadratic filters,” Proc. IEEE, vol. 78, pp. 665-677, Apr. 1990. [I71 G. Ramponi, G. L. Sicuranza, and W. Ukovich, “An optimization approach to the design of nonlinear Volterra filters,” in Proc. EUSIPCO86, Signal Processing III: Theories Applications,, The Hague, The Netherlands, Sept. 1986, pp. 151-154.

963

Sanjit K. Mitra (S’59-M’63-SM’69-F‘74) received the B.Sc. degree (Hons.) degree in physics in 1953 from Uktal University, Cuttack, India, the M.Sc. (Tech.) degree in radio physics and electronics in 1956 from Calcutta University, Calcutta, India, and the M.S. and Ph.D. degrees in electrical engineering from the University of Califomia, Berkeley, in 1960 and 1962, respectively. From June 1962 to June 1965 he was with Comell University, Ithaca, NY, as a Professor of Electrical Engineering. ]He was with AT&T Bell Laboratories, Holmdel, NJ, from June 1965 to Janu.my 1967. He has been with the faculty of the University of California since then, first at the Davis campus and more recently at the Santa Barbara campus as a Professor of Electrical and Computer Engineering, where he served as Chairman of the Department from July 1979 to June 1982. Dr. Mitra served as President of the IEEE Circuits and Systems Society in 1986. He is currently a member of the editorial boards of the Intemational Joumal on Circuits and Systems and Siganl Pvocessing, the Intemational Joumal on Multidimensional Systems and Signal Processing, Signal Processing, and the Joumal of the Franklin Institute. He received the 1973 F. E. Terman Award and the 1985 AT&T Foundation Award of the American Society of Engineering Education, the Education Award of the IEEE Circuits and Systems Society in 1989, and the Distinguished Senior U.S. Scientist Award from the Alexander von Humboldt Foundation of West Germany in 1989. In May 1987, he received an Honorary Doctorate of Technology degree from Tampere University, Tampere, Finland. He is a Fellow of the AAAS and SPIE and is a member of EURASIP and ASEE.

Stefan Thurnhofer (M’96) received the Dip1 -1ng degree in electncal engineering from the University of Erlangen-Nurnberg, Germany, in 1991 and the Ph.D. degree in electrical and computer engineering in 1995 from the University of Califomia, Santa Barbara (UCSB). His research interests at UCSB were in digital signal and image processing, the design and application of nonlinear enhancement methods for digital images using Volterra filters, image interpolation, and error diffusion halftoniug In 1995, he joined Lucent Technologies (formerly, AT&T Bell Laboratone\), Allentown, PA, where he is currently involved in the design of microcontrollers and digital signal processors

Authorized licensed use limited to: INDIAN INSTITUTE OF TECHNOLOGY KHARAGPUR. Downloaded on March 8, 2009 at 22:01 from IEEE Xplore. Restrictions apply.