- Email: [email protected]
Adaptive nonlinear filters for 2D and 3D image enhancement
Signal Processing 67 (1998) 237—254
Adaptive nonlinear filters for 2D and 3D image enhancement Se´bastien Guillon!, Pierre Baylou!, Mohamed Najim!,*, Naamen Keskes" ! Equipe Signal and Image ENSERB and GDR ISIS — CNRS, BP 99, 33 402 Talence Cedex, France " De& partement Image CSTJF — ELF Aquitaine, 64 018 Pau Cedex, France Received 25 June 1997; received in revised form 2 February 1998
Abstract Unsharp masking method is a popular approach for image enhancement, in which a highpass version of an image is added to the original one. This method is easy to run, but is very sensitive to noise. Suppressing noise is generally performed with lowpass filters, and leads to edge blurring. So, an approach which is a combination of a nonlinear lowpass and highpass filters is proposed. These filters are based on an adaptive filter mask. We demonstrate that this approach performs noise reduction as well as edge enhancement. It also improves the contrast enhancement in comparison with other methods. These results are illustrated by processing blurred and noisy images. The method is then extended for 3D data processing and used on 3D seismic images. ( 1998 Elsevier Science B.V. All rights reserved. Zusammenfassung Die Methode zur Unscha¨rfeverdeckung ist ein beliebter Ansatz zur Bildverbesserung, bei dem eine Hochpa{version des Bildes dem Original hinzuaddiert wird. Die Methode ist einfach zu handhaben, aber sehr empfindlich auf Rauschen. Die Rauschunterdru¨ckung wird im allgemeinen mit Tiefpa{filtern durchgefu¨hrt und fu¨hrt zur Kantenverschmierung. Somit wird ein Ansatz vorgeschlagen, der in einer Kombination eines nichtlinearen Tiefpa{- und Hochpa{filters besteht. Diese Filter beruhen auf einer adaptiven Filtermaske. Wir zeigen, da{ dieser Ansatz sowohl eine Rauschreduktion als auch eine Kantenverbesserung vollbringt. Er verbessert im Vergleich zu anderen Methoden auch die Kontrastversta¨rkung. Diese Ergebnisse werden anhand der Verarbeitung unscharfer und verrauschter Bilder illustriert. Die Methode wird sodann auf die Verarbeitung von 3D-Daten erweitert und auf seismische 3D-Bilder angewandt. ( 1998 Elsevier Science B.V. All rights reserved. Re´sume´ La me´thode Unsharp Masking, tre`s utilise´e en re´haussement de contraste, consiste a` ajouter a` l’image une version passe-haut d’elle-meˆme. Cependant cette me´thode, facile a` mettre en oeuvre, est tre`s sensible au bruit. Par ailleurs, la re´duction de bruit est ge´ne´ralement re´alise´e au moyen de filtres passe-bas, ce qui engendre un flou au niveau des contours. Nous proposons dans cet article de combiner deux filtres adaptatifs non-line´aires passe-bas et passe-haut afin de re´aliser dans le meˆme traitement un lissage et un re´haussement de contraste. Nous mettons en e´vidence la robustesse de la me´thode en traitant des images floues et bruite´es. Enfin le filtre est e´tendu au cas tridimensionnel et est applique´ au traitement de blocs sismiques 3D. ( 1998 Elsevier Science B.V. All rights reserved. Keywords: Image enhancement; Edge sharpening; Adaptive filtering; Nonlinear filtering
* Corresponding author. Tel.: 33 56 84 66 74; fax: 33 56 84 84 06; e-mail: [email protected]. 0165-1684/98/$19.00 ( 1998 Elsevier Science B.V. All rights reserved. PII: S 0 1 6 5 - 1 6 8 4 ( 9 8 ) 0 0 0 4 2 - 5
238
S. Guillon et al. / Signal Processing 67 (1998) 237–254
1. Introduction The purpose of image enhancement techniques is to remove noise and to sharpen details in order to improve the visual appearance of an image. A large variety of enhancement filters have been proposed — linear filters, homomorphic filters, L-filters, and the like. In this paper we propose to extend the classical unsharp masking (UM) method [6,7] by using nonlinear filters. The UM method operates by adding a fraction of the highpass filtered version of the input image to the original one (see Fig. 1). This operator is sensitive to noise due to the presence of the linear highpass filter which cannot discriminate signal from noise. Moreover it perceptually enhances image more in dark areas than in lighter ones. Various schemes have therefore been proposed in order to improve the performances of the UM filter. In [8,10—12,17,18], the authors propose the use of quadratic filters instead of a simple linear highpass filter. These filters perform detail and edge enhancement in accordance with the characteristics of human vision. Indeed, enhancement is greater in bright areas than in dark ones in order to respect Weber’s law [6] which states that the just noticeable brightness difference is proportional to the average background brightness. The perceived noise then decreases. In [9,13] Ramponi proposes to use a simple cubic operator which only performs a sharpening action if the processing mask is located across the edge of an object. Based on the same idea, DeFigueiredo has proposed exponential Volterra filters capable of sharpening edges in noisy images [1]. All these techniques present a drawback in having a fixed filter mask for the entire image, which can be
Fig. 1. Unsharp masking (UM) method.
considered as a nonstationary process. In [2,14,16,19,20], the authors propose to use a filter mask that varies according to the local image gradient. The sharpening effect of the edges is achieved, but the method requires a high computation cost due to the large number of iterations needed. In using simultaneously a quadratic filter and a varying filter mask, we propose in [3,5] a family of adaptive quadratic filters extending the unsharp masking operator. The operator formed by a parallel architecture composed of an adaptive linear smoothing and an adaptive quadratic sharpening components, performs noise reduction and edge enhancement. The use of quadratic filters leads to a stronger enhancement of bright pixels than dark ones. Although this treatment is relevant for human perception, it could be detrimental for a post-processing. For instance, in seismic images black and white pixels represent positive and negative reflections, and it is important to enhance them similarly. In this case we propose to simply use a linear adaptive UM method. The proposed approach is developed in two steps. We first define the filter mask as a set of coefficients which quantifies the contribution of each pixel of the support in order to compute the output value. Depending on the type of images we want to process (natural, periodic or pseudo-periodic texture, seismic, etc.), we propose two different ways to compute these coefficients: locally adapted and recursively adapted. The second step consists in combining the pixel values according to this set of coefficients. The aim of this paper is to provide a robust enhancement filter in comparison with other 2D methods, and a 3D extension of this filter in order to process seismic images. The paper is organized as follows. In Section 2 the derivation of the best filter mask according to a performance criterion is presented. In Section 3, the contrast enhancement filter as a combination of a lowpass and a highpass filter is described. The two lowpass and highpass components are then developed and their properties discussed. In Section 4, results of image enhancement are presented in comparison to other approaches. Section 5 is devoted to the 3D extension of the filters, and its
S. Guillon et al. / Signal Processing 67 (1998) 237—254
239
application to the restoration of seismic images. Concluding remarks can be obtained in Section 6.
according to the minimization of a criterion defined later.
2. Filter mask
2.1. Local adaptive mask
Let us define the support S as the set of pixels likely to contribute in the filtering of the current pixel (x,y). Without loss of generality, we consider that S is a rectangular window of size K ]K : 1 2 K !1 K !1 S" (i, j)3 ! 1 ; 1 2 2
The simplest way to define the similarity between two pixels is to compare their grey levels. The approach is based on the observation that variations of gray levels inside a nontextured region are smaller than those between two different regions. Then, for a given location (i, j) of S, we estimate a continuity value D f (x#i, y#j)!f (x, y)D, and the coefficients m are computed by i,j m "C(D f (x#i,y#j)!f (x,y)D), (3) i,j where C is a positive decreasing function, such that C(d)P0 as d increases. A large number of such functions have been proposed in the bibliography for edge-preserving smoothing techniques [2]: inverse gradient filter [20], improved inverse gradient filter [19] and adaptive Gaussian weighted filter [16]. Hence a possible choice for C is
G C
C
D
K !1 K !1 ] ! 2 ; 2 2 2
DH
.
(1)
When scanning the image, the pixels included in the mask may belong to different regions. The principle of adaptive filtering then consists in adjusting the contribution of each pixel of S according to their similarity with the current pixel. For this purpose, we define the filter mask M as a set of coefficients m (in the range [0,1]) associated with i,j the support S: M"Mm 3[0,1]D(i, j)3SN. (2) i,j Each coefficient m is viewed as a level of confidence i,j for the pixel (x#i, y#j) to belong to the same region as the current pixel (x,y), i.e. m P1 for i,j pixels similar to the current pixel, respectively m P0 for the others. i,j In order to compute the set of coefficients M, a criterion must be defined to estimate the similarity between two pixels. A first criterion could be defined in comparing the two corresponding grey levels. The deviation between the two pixel values is calculated, and a value between 0 and 1 is assigned to m such that i,j m "1 when the deviation is zero, and m "0 for i,j i,j large deviation. This method is referred to a local adaptive approach. This estimation is fairly appropriate for natural images, for which each pixel to be filtered needs a specific mask. A second approach, devoted to periodic or pseudo-periodic textured images, is also presented. For this kind of images, the filter mask is stable or slow spatially-varying. Then we propose a recursive approach, in which the mask M evolves iteratively
A
B
( f (x#i,y#j)!f (x,y))2 m "exp ! , i,j 2d2
(4)
where d is a parameter controlling the width of the curve or the effective degree of similarity between the two positions. Its choice will be discussed in Section 4. The function has proved, through a large number of simulations, to give better results than those presented in [19,20]. Intuitively, from Eq. (4) we see that pixels belonging to the same region as the pixel (x,y) will have larger weighting coefficients than those located outside the region. Remark. In order to reduce the computation cost due to the Gaussian term in the above equations, note that the possible continuity values are limited to 255 values, for an 8-bit coded image. In practice a table can be used to estimate the m values. i,j 2.2. Recursive adaptive mask The idea of adapting recursively the filter mask is due to Salembier who derived an adaptive approach
240
S. Guillon et al. / Signal Processing 67 (1998) 237–254
for rank order-based filters [15] to process periodic or quasi-periodic textured images. We propose here another adaptation process of this method. The method, proposed by Salembier to define an adapted unweighted filter mask, consists in defining a search area, and assigning a weighting coefficient to each possible location. The current filter mask support is obtained by thresholding the set of coefficients — if a coefficient is greater than the threshold, the corresponding location is considered as belonging to the unweighted filter mask. An adaptation process is used to optimize (according to a criterion) the filter mask support by modifying the set of continuous values of the search area. In the approach we propose, a weighted filter mask is used instead of a binary one. All the pixels (x#i, y#j) belonging to a predefined support S, are taken into account according to their weight m — no arbitrary thresholding is necessary to i,j obtain a binary mask. The algorithm then operates in two steps: it first computes the coefficients m , i,j and secondly implements the linear filtering. The structure of the filter is given in Fig. 2 where the coefficients are adaptively modified by minimizing a prediction error. The predicted value p(x,y) we propose is a weighted mean over S*"SCM(0,0)N: f (x#i,y#j) + *m p(x,y)" (i,j)|S i,j . (5) m + (i,j)|S* i,j The optimality criterion we consider in this paper is the mean square error (MSE): J"E[(p(x,y)!f (x,y))2],
(6)
where E[)] is the mathematical expectation. Similar results are obtained with a mean absolute error (MAE) criterion. The minimization of the prediction error is performed with a steepest descent algorithm [21], chosen for the simplicity of the LMS algorithm, and its ability to deal with nonstationary signals. The minimization process is then reached by computing an estimate of the gradient of the criterion J with respect to the set of parameters M. If (mk ) i,j i,j|S is the current set of coefficients at step k, the next estimate at step k#1 is computed with the
Fig. 2. Filter structure.
following formula:
A
B
LJ mk`1"W mk !j , (7) i,j i,j Lm i,j where j is a parameter which controls the convergence of the algorithm, and W is a function maintaining the m values in the range [0,1] in order to i,j satisfy the definition of the filter mask M (Eq. (2)). The W function can be a simple clamping or an affine correction as follows: k !min k , mk`1" i,j i,j max !min k k where LJ k "mk !j , i,j i,j Lm i,j
(8)
(9)
min " min Mk N, k i,j (i,j)|S
(10)
max "max Mk N, (11) k i,j (i,j)|S In the derivation of the LMS algorithm, the gradient estimate of the criterion involved in Eq. (6) is computed by simply approximating the expectation with the instantaneous value. This leads to the following estimate: LJ LE[(p(x,y)!f (x,y))2] " Lm Lm i,j i,j Lp(x,y) "2E (p(x,y)!f (x,y)) Lm i,j Lp(x,y) K2(p(x,y)!f (x,y)) , Lm i,j
C
D (12)
S. Guillon et al. / Signal Processing 67 (1998) 237—254
241
with m !+ m f (x#r, y#s) Lp(x,y) f (x#i, y#j)+ p(x,y)!f (x#i, y#j) (r,s)|S* r,s (r,s)|S* r,s " "! . m )2 m (+ Lm + (r,s)|S* r,s i,j (r,s)|S* r,s So we obtain the following updating algorithm:
C
mk`1"t mk #2j(p(x,y)!f (x,y)) i,j i,j
A
BD
p(x,y)!f (x#i, y#j) ] . (14) mk + (r,s)|S* r,s As Eq. (14) is a nonlinear function of the coefficients m , we cannot easily carry out its convergence i,j properties. However, two well-known properties of the LMS algorithm [21] are still valid in this case: f the coefficients m exhibit an exponential coni,j vergence, with a time constant q inversely proportional to the parameter j: 1 qJ ; j f
(15)
the final misadjustment of the coefficients increases with j, and the mean square estimation error e is proportional to j: mi,j (16) E[e2 ]Jj. mi,j So the choice of j results from a compromise between the convergence speed (necessary in the case of non-stationary images) and the misadjustment of the coefficients m . i,j One can observe in the updating equation (14) that as p(x,y) is a predicted value of f (x,y), the set M evolves in order to improve the estimation p(x,y) of f (x,y). Indeed f if p(x,y) underestimates f (x,y) (i.e. p(x,y)( f (x,y)), then the values m for which f (x#i,y#j)' i,j p(x,y) are increased; f if p(x,y) overestimates f (x,y) (i.e. p(x,y)' f (x,y)), then the values m for which f (x#i,y#j)( i,j p(x,y) are increased. By copying this procedure, one can derive several updating equations which increase coefficients of pixels similar to the current pixel, and decrease other coefficients. For example, by introducing a similarity index r (in the range [!1,1]) between i,j the two pixels (x,y) and (x#i,y#j), which is max-
(13)
imum for two similar values and minimum for two different values, we propose the general following updating equation: (17) mk`1"W[mk #jr ]. i,j i,j i,j In the last section we will provide an example of such an updating function for the restoration of seismic images.
3. Contrast enhancement filter In order to improve the performances of the UM method, we propose to combine a lowpass and a highpass filter (see Fig. 3). The lowpass component g is used to smooth data in homogeneous areas, L while the highpass g is chosen as an edge detector H in order to sharpen details in the image. The filter output g(x,y) is given by g(x,y)"g (x,y)#ag (x,y), L H where
(18)
g (x,y)" + h f (x#i,y#j), L i,j (i,j)|S
(19)
g (x,y)" + w f (x#i,y#j). (20) H i,j (i,j)|S The coefficient a is a factor driving the edge enhancement effect. The operator we propose is then defined by the choice of the coefficients h and i,j w which are different combinations of the m i,j i,j
Fig. 3. Enhancement filter structure.
242
S. Guillon et al. / Signal Processing 67 (1998) 237–254
presented in the previous section. As reported in [10], to preserve the expected output of a uniform luminance input, the following conditions must hold:
S of pixels the most similar to the current pixel. 1 This threshold could be set arbitrary to a constant value (for example ¹"0.5), or be estimated adaptively by taking
+ h "1 and + w "0. (21) i,j i,j (i,j)|S (i,j)|S In the following we first give the expression of the lowpass and highpass filters, and finally propose an adaptive estimation of the scaling factor a.
1 ¹"mN " + m . (26) i,j Card(S) (i,j)|S With this value the class S is the set of pixels 1 having the highest levels of confidence.
3.1. Lowpass filter
3.2. Highpass filter: ¸aplacian like enhancement (¸¸E)
As the purpose of lowpass filtering is to reduce noise in the image without smoothing details or edges, we consider the following lowpass filter g (x,y): L g (x,y)" + h f (x#i, y#j ), L i,j (i,j)|S m i,j with h " . (22) i,j + m (k,l)|S k,l The condition (21) can be easily checked: m i,j "1. (23) + h " + i,j + m (i,j)|S (i,j)|S (k,l)|S k,l It results from the definition of the mask M, that g (x,y) is a weighted mean over the set of pixels the L most similar to (x,y) (pixels having high m values). i,j So the lowpass filter effectively performs noise reduction without edge smoothing. This filter belongs to the class of edge preserving smoothing filters (EPSF) widely studied in [2,14,16,19,20]. The performances of the lowpass filter can be improved by defining the following submask S : 1 S "M(i, j)3SDm '¹N, (24) 1 i,j where ¹3[0,1] is a threshold, and by performing the new lowpass filter: g1(x, y)" + h1 f (x#i,y#j), L i,j (i,j)|S1 m i,j with h1 " . (25) i,j + m (k,l)|S1 k,l According to the choice of the threshold ¹, g1(x,y) L becomes a weighted mean restricted over a class
In this section, we propose a new highpass filter defined by a linear combination of the coefficients of the filter mask M. We show that it is equivalent to a Laplacian filter, but that it exhibits an improved noise robustness. The highpass filter g (x,y) we propose is defined H by g (x,y)" + w f (x#i,y#j), H i,j (i,j)|S w "m !mN , i,j i,j with mN " 1 + m . C!3$(S) i,j (i,j)|S The set of coefficients W"(w ) verifies Eq. i,j + w " + (m !mN )"0. i,j i,j (i,j)|S (i,j)|S To point out the behaviour of the highpass ponent, we write g (x,y) as follows: H
G
(27) (21): (28) com-
g (x,y)" + (m !mN ) f (x#i,y#j) H i,j (i,j)|S " + (m !mN ) f (x#i,y#j) i,j (i,j)|S1 # + (m !mN ) f (x#i,y#j) i,j (i,j)|S2 " + Dm !mN D f (x#i,y#j) i,j (i,j)|S1 ! + Dm !mN D f (x#i,y#j), i,j (i,j)|S2
(29)
S. Guillon et al. / Signal Processing 67 (1998) 237—254
A
n(N!n) + f (x#i, y#j) (i,j)|S1 K N n
with S "M(i, j)3SDm *mN N 1 i,j set of pixels similar to current pixel,
B
f (x#i, y#j) + ! (i,j)|S2 N!n
S "M(i, j)3SDm (mN N 2 i,j set of pixels different from current pixel. So g (x,y) can be interpreted as the local difference H of pixel values belonging to two different regions — the ones belonging to S which are similar to the 1 current pixel, and the others belonging to S which 2 are different from the current pixel. The filter response can be studied in the following cases: f homogeneous area If the current pixel (x,y) is located on a quasi uniform area, then we can state that f (x#i,y#j)+f (x,y)Nm K1 ∀(i, j)3S i,j NmN K1
f
NDm !mN DK0 ∀(i, j)3S i,j Ng (x,y)+0. (30) H The output of the highpass filter is then negligible. presence of an edge If the pixel (x,y) is located across an edge (Fig. 4), the two classes S and S are characterized by 1 2 (i, j)3S Nm K1, 1 i,j (31) (i, j)3S Nm K0. 2 i,j So considering that Card(S)"N, Card(S )"n 1 and Card(S )"N!n, we have mN Kn/N and 2 n g (x,y)K + 1! f (x#i,y#j) H N (i,j)|S1 n ! + 0! f (x#i,y#j) N (i,j)|S2
A
A
B
B
243
n(N!n) K ( fM !fM ), @S1 @S2 N
(32)
with fM and fM the mean values of f estimated @S2 @S1 over S and S , respectively. Thus we easily show 1 2 that if (x,y) is located at the bottom of the edge, then g (x,y)(0 as fM (fM . If (x,y) is located at H @S1 @S2 the top of the edge, then g (x,y)'0 as fM 'fM . H @S1 @S2 And finally if (x,y) is located along the edge then g (x,y)"0 as fM "fM . H @S1 @S2 Globally the highpass filter is then equivalent to a classical Laplacian operator, but with much less noise sensitivity as we will see further (Section 4). More generally, g (x,y) is equivalent to a direcH tional second order derivative considered orthogonally to the edge. For example, a support S of size 3]3 located across an edge gives the following filter mask:
(33) with a a parameter which depends on the way the filter mask is implemented.
Fig. 4. Laplacian behavior.
244
S. Guillon et al. / Signal Processing 67 (1998) 237–254
Fig. 5. Filter structure with an adaptive a. Fig. 6. Adaptive a.
3.3. Scaling factor a As introduced in Eq. (18), a is a factor driving the enhancement effect. Choosing a constant value for a has the drawback of enhancing sharp edges more than smoothed ones because the enhancement is proportional to the sharpness of the detected edge. This implies an overshooting effect on the image (false black or white lines along edges), which could result in an unpleasant enhancement effect. This can be solved by applying a local adaptation of a depending on the sharpness of the edge. As the output of the highpass filter measures this sharpness, we use the value g (x,y) to control H a (Fig. 5): a"/(Dg (x,y)D). H
(34)
The enhancement function /()) must be continuous (/ : R`PR`) and satisfy f lim /(x)"0: a small g value is considered x?0 H as a false detection and no enhancement is required. f lim /(x)"0: a high g value corresponds x?= H to a sharp edge, and no enhancement is required. Lagrangian functions (Eq. (35)) satisfy all these constraints: a"/(Dg (x,y)D)"CDg (x,y)Dexp(!cDg (x,y)D). H H H
(35)
We can also use a simple piecewise linear approximation of this function (Fig. 6). Threshold ¹ controls noise removal and thre1 shold ¹ controls the overshooting effect. The 2 simulations presented in the last section are carried out with ¹ "5 and ¹ "50 (for 8-bit coded 1 2 images).
3.4. Iterative application of the enhancement filter The new proposed filter performs simultaneously edge enhancement and region smoothing. However, according to the level of noise and blur in the image, a single processing of the image could be not sufficient to achieve enough enhancement. We suggest in this case to apply iteratively the enhancement filter. Let g(x,y,n) be the filtered image at iteration n, then the filtering process becomes g(x,y,0)"f (x,y),
(36) g(x,y,n#1)"g (x,y,n)#ag (x,y,n) ∀n*0. L H This process is equivalent to that proposed in [14], but needs a smaller number of iterations — less than five iterations are sufficient in practice. We provide an example in the next section.
4. Application to image enhancement and filtering In this section we will describe some practical examples which demonstrate the performances of these new filters to restore blurred or noisy images. Simulations are performed with different classes of images — synthetic, natural and textured.
4.1. Application of the local filter In this section, the coefficients are calculated from the deviation between the current pixel value and the values of the pixels included in support S (Eq. (4)). The parameter d has to be set at first.
S. Guillon et al. / Signal Processing 67 (1998) 237—254
4.1.1. Choice of d In order to state how to choose the parameter d involved in the filter mask, we first study the behavior of the new filter on a noisy region with or without the presence of an edge. First, if the pixel (x,y) is located on a uniform area (mean value f ) corrupted by a white Gaussian 0 noise v(x,y) of variance p2: f (x,y)"f #v(x,y). (37) 0 Then the set of coefficients of the filter mask is defined by ( f (x,y)!f (x#i,y#j))2 m "exp! i,j 2d2 (v(x,y)!v(x#i,y#j))2 "exp! . 2d2
(38)
Assuming that v(x,y) is a Gaussian random variable, m is a realisation of a random variable m whose i,j density function p (m) can be expressed by m D m(D2~2)@2 d p (m)" m3[0,1] with D" . m p J2n J!ln m (39) Then, the mean mN and the variance p2 of these m coefficients can be obtained:
P
mN "E[m]" and
1
0
D mp (m)dm" m J2#D2
(40)
P
1 m2p (m)dm!mN 2 m 0 D D2 " ! . J4#D2 2#D2
p2 "E[(m!mN )2]" m
(41)
We can state from relations (40) and (41) that if we take D'2, then mN '0.8 and p2 )0.04, which m implies Nm KmN i,j
∀(i, j)3S
+ m ) g(x#i,y#j) Ng (x,y)" (i,j)|S i,j PB + m (i,j)|S i,j 1 K + g(x#i,y#j) Card(S) (i,j)|S
(42)
Ng (x,y)" + (m !mN ) ) g(x#i,y#j)K0. PH i,j (i,j)|S
245
This means that if d verifies d'2p, then the enhancement filter is equivalent to a smoothing filter, and the output variance is given by p2 p2" . g Card(S)
(43)
In the second case in which the pixel (x,y) is located across an edge, the choice of d depends on the height h of this edge. If a and b are the grey levels on the two sides of the edge, and if f (x,y)"a we can state that (i, j)3S Q f (x#i,y#j)Ka N m "m K1, 1 i,j 1 (i, j)3S Q f (x#i,y#j)Kb 2 (b!a)2 N m "m Kexp! i,j 2 2d2 h2 "exp! . (44) 2d2 The enhancement effect of the filter is effective only if the two classes S and S are well separated, i.e. if 1 2 m @m , which implies 2 1 (b!a)2 exp! @1 Q d@h. (45) 2d2 Finally, we can state that the choice of the parameter d must satisfy the following condition: 2p(d(h , (46) .*/ where p2 is the noise variance, and h the smallest .*/ height of edge we want to enhance. 4.1.2. Application In order to illustrate the behavior and the noise robustness of this enhancement filter, we compare the response of the UM and LLE methods. We first consider a blurred edge (h"110), and we compare results obtained by the highpass component: LLE (3]3, d"20) and Laplacian (3]3). In Fig. 7, one can observe that these two components are similar, and finally the edge enhancement is equivalent (Fig. 8). We then degrade the synthetic edge by adding a white Gaussian noise of variance p2"16. Results of the highpass components are provided in Fig. 9 — the support S is of size 5]5, and d"40 for the LLE. These results show the noise robustness of the LLE operator in comparison with the Laplacian.
246
S. Guillon et al. / Signal Processing 67 (1998) 237–254
Fig. 7. Comparison of the Laplacian behavior.
Fig. 8. Comparison of edge enhancement.
S. Guillon et al. / Signal Processing 67 (1998) 237—254
247
Fig. 9. Comparison of the Laplacian behavior on a noisy edge.
As expected, the new enhancement filter allows noise smoothing in uniform areas in the same time than edge enhancement (Fig. 10). In the second application, we propose to enhance the image of boat presented in Fig. 11. If we exercise the quadratic operators proposed by Mitra and Ramponi [8,10] on this image, we observe noise amplification in background area (such as the sky). The cubic operator proposed by Ramponi [9] is more robust against noise, but one can observe that some continuous edges are discarded. In Figs. 12 and 13 we provide the results obtained with the Mitra quadratic operator, the Ramponi cubic filter and the LLE operator, respectively. In these images, the contrast enhancement effect of the LLE operator is easily perceived without exhibiting noise amplification in background areas. An adaptive scaling factor a allows enhancement without overshoot — although the image obtained with an adaptive a is perceptually not as good as that obtained with a constant a, the absence of overshoot allows to use this image for a postprocessing.
4.2. Application of the recursive filter Some examples of enhancement of periodic or quasi-periodic textured images are presented to point out the interest of the recursive filter mask. Indeed, on such images, the optimal filter mask will be constant or slow spatially-varying. An adaptive evaluation of the filter mask M will then be more robust than a local one — if the pixel to be filtered is corrupted by noise, the coefficients of the filter mask will not be affected due to the high memory effect of M. The image to be filtered is the ‘canvas’ image presented in Fig. 14. If we process such an image by a quadratic or a cubic operator (Fig. 15: Ramponi’s cubic filter [9] with d"0.0003, Mitra’s quadratic filter [8] with a"512), we observe a noise amplification which gives unpleasant filtered images. On the contrary, results obtained with a 3]3 LLE filter appear more convincing in terms of contrast enhancement (Fig. 16: LLE method) — the texture is well enhanced without noise amplification. The image is scanned along columns so that the coefficients quickly converge.
248
S. Guillon et al. / Signal Processing 67 (1998) 237–254
Fig. 10. Comparison of edge enhancement on a noisy edge.
blurring. As the LLE operator combines a lowpass filter to remove noise and a highpass filter to restore edges, restoration with this operator will be more robust. We have iteratively filtered the noisy image with a"0.4, j"0.01 and S"[5]5] (Fig. 17). After four iterations the image becomes nearly binary with a good restoration of its main features, whose perception are difficult in the original image.
5. 3D extension and application to seismic images
Fig. 11. ‘Boat’ image.
In the second application, we propose to restore the ‘canvas’ image corrupted by a zero-mean Gaussian white noise with standard deviation p2"60. Classically additive Gaussian noise is removed by applying low pass filters, but at the expense of edge
The main purpose in seismic data analysis is the detection of ‘geological horizons’ separating homogeneous layers of rocks, sediments, etc. A 2D seismic image (Fig. 19) is composed of seismic traces, sinusoid-like waveforms, which are a record of reflected waves arising from impedance contrasts between strata. If a cycle is correlated laterally across many seismic traces on an image, it is called a seismic horizon. Filtering seismic images then consists in restoring these seismic horizons. As a seismic image is often
S. Guillon et al. / Signal Processing 67 (1998) 237—254
249
Fig. 12. Image enhancement.
Fig. 13. Local zoom.
K ]K ]K (corresponding support S). 1 2 3
G
C
D
K !1 K !1 S" (i,j,k)3 ! 1 ; 1 2 2
C C
D DH
K !1 K !1 ] ! 2 ; 2 2 2 K !1 K !1 ] ! 3 ; 3 2 2 M"Mm
,
3[0,1]D(i, j, k)3SN. i,j,k
(47) (48)
The output g(x,y,z) of the filter then becomes Fig. 14. ‘Canvas’ image.
extracted from a 3D bloc of data, the filtering process can be implemented over a 2D or 3D support. In the case of 3D processing, the filter mask M is naturally extended as a window of size
g(x,y,z)"g (x,y,z)#ag (x,y,z), L H
(49)
with the component g and g defined by a simple L H 3D extension: m f (x#i,y#j,z#k) i,j,k g (x,y,z)" + , L + m (i,j,k)|S i,j,k (i,j,k)|S
250
S. Guillon et al. / Signal Processing 67 (1998) 237–254
Fig. 15. Quadratic enhancement.
g (x,y,z)" + (m !mN ) f (x#i,y#j,z#k), H i,j,k (i,j,k)|S 1 with mN " + m . (50) i,j,k Card(S) (i,j,k)|S In the following, we propose local adaptive and recursively adaptive methods, specifically adapted to seismic images, due to the way of computing the set of coefficients M. To simplify the presentation, the equations are given in the 2D case, the 3D case being a simple extension by adding the third coordinate z.
5.1. Seismic filter mask In Section 2, we defined various criteria of similarity between two pixels in order to compute the coefficients m . Here we propose to improve the i,j similarity criterion for seismic images. Indeed two pixels could be said to be similar if they belong to the same seismic horizon, which means that their corresponding seismic traces are highly correlated. A seismic trace is a column of the image, and its representation is a sinusoid-like waveform (Fig. 18: intensity evolution along a column of the image). So for each position (x,y) we define the vector seismic trace by (x,y)"[ f (x!k,y)] . 53 k|*~N@2§N@2+
(51)
Fig. 16. LLE enhancement.
The similarity of two pixels is then locally estimated by measuring the correlation factor o of their two corresponding traces, which gives in the 2D case: 1 1 m " # o( (x,y), (x#i,y#j)) i,j 2 2 53 53 1 1 T (x,y) ) (x#i,y#j) 53 53 " # . 2 2 E (x,y)E ) E (x#i,y#j)E 53 53
(52)
S. Guillon et al. / Signal Processing 67 (1998) 237—254
251
Fig. 17. Iterative enhancement.
A second approach for computing these coefficients can be derived from the recursive scheme exposed in Section 2.2. Indeed we assume that the generic updating equation expressed by mk`1"t[mk #jr ], (53) i,j i,j i,j where r is a measurement of similarity betwen two pixels respecting r "#1, if (x,y) and (x#i,y#j) are similar, i,j r "!1, if (x,y) and (x#i,y#j) (54) i,j are widely different.
Fig. 18. Seismic trace.
Thus we could assume that m 3[0,1], and that i,j m "1 only if (x,y)"b (x#i,y#j) (b'0). As i,j 53 53 it takes into account the vicinity of the pixels, the estimation of the coefficients m is more robust to i,j noise than the one proposed in Section 2.1.
As the correlation factor o (Eq. (52)) of two seismic traces satisfies these conditions, we propose to define the updating equation by
A
B
T (x,y) ) (x#i,y#j) 53 mk`1"t mk #j 53 . i,j i,j E (x,y)E ) E (x#i,y#j)E 53 53
(55)
Thus we obtain an adaptive estimation of the set of coefficients M which tends to increase the weights
252
S. Guillon et al. / Signal Processing 67 (1998) 237–254
of the positions the most correlated to the current position. The parameter j drives the memory of the filter mask.
5.2. Application The two proposed schemes we have presented are used to process the seismic image shown in Fig. 19. We first process this image with a 3D locally adapted mask, whose weights are evaluated with the correlation factor of the seismic traces (Eq. (52)), S"[5]5]5], a seismic trace of size 11 (N"5) and an enhancement factor a"0.1. The result is presented in Fig. 20. To evaluate the filtering effect of our approach on this class of images, we should focus our attention on the main features of these images: seismic hor-
izons, channels (located in the upper left of the image) and faults (right side of the image). Important features for expert interpretation are preserved while an efficient cancellation of noise is performed. We compare this result with the filtering using a recursive mask, where the coefficients evolve according to the correlation factor of the seismic traces (Eq. (55)). We scan the image in the column direction so that the filter mask is slow spatiallyvarying. The image is processed with an S"[5]5]5], a vector seismic trace of size 11, an updating coefficient j"0.01, and an enhancement factor a"0.3. The result is given in Fig. 21. The evaluation of the result is still done by focusing our attention to the features of this image: seismic horizons appear more regular than those resulting from the local filter (Fig. 20), but the channel and the fault are not as well preserved as with the local filter. Finally, the two approaches (local and recursive) can be combined if we are able to detect the position of the two main features (fault and channel). If we note d(x,y) the probability function of such a detection (d(x,y)P1 if a feature is detected), then the output image is obtained by
Fig. 19. Seismic image.
g(x,y)"g (x,y) ) d(x,y)#g (x,y) ) (1!d(x,y)), (56) -0# 3%# with g and g respectively being the local and -0# 3%# recursive filtering. From Fig. 22, one can observe that the main features are well preserved, and horizons regularized. This method combines the advantages of the two previous ones.
Fig. 20. Seismic filtered image with a local mask.
Fig. 21. Seismic filtered image with a recursive mask.
S. Guillon et al. / Signal Processing 67 (1998) 237—254
253
Notations f ()) g()) (x,y) (x,y,z) S (i,j) (i,j,k) Fig. 22. Seismic filtered image combining the local and recursive mask.
M"(m ) (or (m )) i,j i,j,k ()M ) E[)]
The detection d(x,y) we used is based on a measure of dispersion of local gradients [4] not described in this paper.
6. Conclusion In this paper new nonlinear image enhancement filters have been presented. They are based on the choice of a locally-adapted filter mask or a recursive adaptive process for computing the mask coefficients. The filter mask defines the set of surrounding pixels most closely correlated to the current pixel. Then two new enhancement techniques, called Laplacian like enhancement (LLE1 and LLE2), have been introduced as an improvement of the classical unsharp masking (UM) method. We have shown, through tests on synthetic and real images, the improvements brought by these algorithms: noise reduction and edge enhancement are performed simultaneously. Comparisons with other techniques have also been provided. Finally, these new approaches, including a 3D extension, have been applied to practical applications dealing with seismic images restoration. From a computational point of view, the use of a locally-adapted filter mask is not more costly than other UM algorithms (parallel implementation is possible), while the adaptive process for computing the coefficients leads to a more costly algorithm.
input image output image position of a pixel on an image position of a voxel on a 3D bloc of data filtering support position of a pixel on a 2D support position of a voxel on a 3D support set of levels of confidence (respectively 2D or 3D) mean operator mathematical expectation
Acknowledgements The authors wish to thank ELF Aquitaine for having provided a large number of seismic images, and for partially supporting this work. They are also indebted to the academic partners of the joint project CNRS — ELF Aquitaine. They also thank the anonymous reviewers for their helpful comments that aided a lot in making this paper more clear and complete.
References [1] R.J.P. de Figueiredo, S.C. Matz, Exponential nonlinear Volterra filters for contrast sharpening in noisy images, in: Proc. ICASSP’96, Vol. 4, pp. 2263—2266. [2] J.M.H. Du Buf, T.G. Campbell, A quantitative comparison of edge-preserving smoothing techniques, Signal Processing 21 (4) (1990) 289—301. [3] S. Guillon, P. Baylou, New noncausal adaptive quadratic filters for image filtering and contrast enhancement, in: Proc. EUSIPCO’96, Trieste, Italy, Vol. 1, 1996, pp. 583—586. [4] S. Guillon, P. Baylou, N. Keskes, Detection and/or determination method of features related to singular points, Submitted for a French patent No. 9708602, INPI (Institut National de la Protection Industrielle), July 1997. [5] S. Guillon, P. Baylou, M. Najim, Robust nonlinear contrast enhancement filters, in: Proc. ICIP’96, Lausanne, Switzerland, Vol. 1, 1996, pp. 757—760.
254
S. Guillon et al. / Signal Processing 67 (1998) 237–254
[6] A.K. Jain, Fundamentals of Digital Image Processing, Prentice-Hall, Englewood Cliffs, NJ, 1989. [7] J.S. Lee, Digital image enhancement and noise filtering by use of local statistics, IEEE Trans. Pattern Anal. Machine Intell. 2 (1980) 165—168. [8] S.K. Mitra, H. Li, I. Lin, T. Yu, A new class of nonlinear filters for image enhancement, Proc. ICASSP’91, pp. 2525—2528. [9] G. Ramponi, A simple cubic operator for sharpening an image, Proc. IEEE Workshop on NSIP, Thessaloniki, Greece, 20—22 June 1995, pp. 963—966. [10] G. Ramponi, P. Fontanot, Enhancing document images with a quadratic filter, Signal Processing 33 (1) (1993) 23—24. [11] G. Ramponi, G.L. Sicuranza, Quadratic digital filters for image processing, IEEE Trans. ASSP 36 (6) (June 1988) 1263—1285. [12] G. Ramponi, G.L. Sicuranza, Image sharpening using a polynomial filter, Proc. ECCTD’93, Davos, Switzerland, pp. 1431—1436. [13] G. Ramponi, N. Strobel, S. Mitra, T.H. Yu, Nonlinear unsharp masking methods for image contrast enhancement, Journal of Electronic Imaging 5 (3) (July 1996) 353—366.
[14] P. Saint-Marc, J.S. Chen, G. Medioni, Adaptive smoothing: a general tool for early vision, IEEE Trans. Pattern Anal. Machine Intell. 13 (6) (June 1991) 514—529. [15] P. Salembier, Adaptive rank order based filters, Signal Processing 27 (1) (1992) 1—25. [16] M. Spann, A. Nieminen, Adaptive gaussian weighted filtering for image segmentation, Pattern Recognition Letters 8 (1988) 251—255. [17] S. Thurnhofer, Quadratic Volterra filters for edge enhancement and their application in image processing, Ph.D. Thesis, University of California, Santa Barbara (1994). [18] A. Vanzo, G. Ramponi, G.L. Sicuranza, An image enhancement technique using polynomial filters, in: Proc. ICIP’94, Austin, TX, 13—16 November 1994. [19] D. Wang, O. Wang, A weighted averaging method for image smoothing, Proc. IEEE Conf. Com. Vision and Pattern Recognition, Miami Beach, Florida, USA, 22—26 June 1986, pp. 981—983. [20] D. Wang, A. Vagnucci, C.C. Li, Gradient inverse weighted smoothing scheme and the evaluation of its performances, CGIP 15 (1981) 167—181. [21] B. Widrow, S.D. Stearns, Adaptive Signal Processing, Prentice-Hall, Englewood Cliffs, NJ, 1985.
Adaptive nonlinear filters for 2D and 3D image enhancement Se´bastien Guillon!, Pierre Baylou!, Mohamed Najim!,*, Naamen Keskes" ! Equipe Signal and Image ENSERB and GDR ISIS — CNRS, BP 99, 33 402 Talence Cedex, France " De& partement Image CSTJF — ELF Aquitaine, 64 018 Pau Cedex, France Received 25 June 1997; received in revised form 2 February 1998
Abstract Unsharp masking method is a popular approach for image enhancement, in which a highpass version of an image is added to the original one. This method is easy to run, but is very sensitive to noise. Suppressing noise is generally performed with lowpass filters, and leads to edge blurring. So, an approach which is a combination of a nonlinear lowpass and highpass filters is proposed. These filters are based on an adaptive filter mask. We demonstrate that this approach performs noise reduction as well as edge enhancement. It also improves the contrast enhancement in comparison with other methods. These results are illustrated by processing blurred and noisy images. The method is then extended for 3D data processing and used on 3D seismic images. ( 1998 Elsevier Science B.V. All rights reserved. Zusammenfassung Die Methode zur Unscha¨rfeverdeckung ist ein beliebter Ansatz zur Bildverbesserung, bei dem eine Hochpa{version des Bildes dem Original hinzuaddiert wird. Die Methode ist einfach zu handhaben, aber sehr empfindlich auf Rauschen. Die Rauschunterdru¨ckung wird im allgemeinen mit Tiefpa{filtern durchgefu¨hrt und fu¨hrt zur Kantenverschmierung. Somit wird ein Ansatz vorgeschlagen, der in einer Kombination eines nichtlinearen Tiefpa{- und Hochpa{filters besteht. Diese Filter beruhen auf einer adaptiven Filtermaske. Wir zeigen, da{ dieser Ansatz sowohl eine Rauschreduktion als auch eine Kantenverbesserung vollbringt. Er verbessert im Vergleich zu anderen Methoden auch die Kontrastversta¨rkung. Diese Ergebnisse werden anhand der Verarbeitung unscharfer und verrauschter Bilder illustriert. Die Methode wird sodann auf die Verarbeitung von 3D-Daten erweitert und auf seismische 3D-Bilder angewandt. ( 1998 Elsevier Science B.V. All rights reserved. Re´sume´ La me´thode Unsharp Masking, tre`s utilise´e en re´haussement de contraste, consiste a` ajouter a` l’image une version passe-haut d’elle-meˆme. Cependant cette me´thode, facile a` mettre en oeuvre, est tre`s sensible au bruit. Par ailleurs, la re´duction de bruit est ge´ne´ralement re´alise´e au moyen de filtres passe-bas, ce qui engendre un flou au niveau des contours. Nous proposons dans cet article de combiner deux filtres adaptatifs non-line´aires passe-bas et passe-haut afin de re´aliser dans le meˆme traitement un lissage et un re´haussement de contraste. Nous mettons en e´vidence la robustesse de la me´thode en traitant des images floues et bruite´es. Enfin le filtre est e´tendu au cas tridimensionnel et est applique´ au traitement de blocs sismiques 3D. ( 1998 Elsevier Science B.V. All rights reserved. Keywords: Image enhancement; Edge sharpening; Adaptive filtering; Nonlinear filtering
* Corresponding author. Tel.: 33 56 84 66 74; fax: 33 56 84 84 06; e-mail: [email protected]. 0165-1684/98/$19.00 ( 1998 Elsevier Science B.V. All rights reserved. PII: S 0 1 6 5 - 1 6 8 4 ( 9 8 ) 0 0 0 4 2 - 5
238
S. Guillon et al. / Signal Processing 67 (1998) 237–254
1. Introduction The purpose of image enhancement techniques is to remove noise and to sharpen details in order to improve the visual appearance of an image. A large variety of enhancement filters have been proposed — linear filters, homomorphic filters, L-filters, and the like. In this paper we propose to extend the classical unsharp masking (UM) method [6,7] by using nonlinear filters. The UM method operates by adding a fraction of the highpass filtered version of the input image to the original one (see Fig. 1). This operator is sensitive to noise due to the presence of the linear highpass filter which cannot discriminate signal from noise. Moreover it perceptually enhances image more in dark areas than in lighter ones. Various schemes have therefore been proposed in order to improve the performances of the UM filter. In [8,10—12,17,18], the authors propose the use of quadratic filters instead of a simple linear highpass filter. These filters perform detail and edge enhancement in accordance with the characteristics of human vision. Indeed, enhancement is greater in bright areas than in dark ones in order to respect Weber’s law [6] which states that the just noticeable brightness difference is proportional to the average background brightness. The perceived noise then decreases. In [9,13] Ramponi proposes to use a simple cubic operator which only performs a sharpening action if the processing mask is located across the edge of an object. Based on the same idea, DeFigueiredo has proposed exponential Volterra filters capable of sharpening edges in noisy images [1]. All these techniques present a drawback in having a fixed filter mask for the entire image, which can be
Fig. 1. Unsharp masking (UM) method.
considered as a nonstationary process. In [2,14,16,19,20], the authors propose to use a filter mask that varies according to the local image gradient. The sharpening effect of the edges is achieved, but the method requires a high computation cost due to the large number of iterations needed. In using simultaneously a quadratic filter and a varying filter mask, we propose in [3,5] a family of adaptive quadratic filters extending the unsharp masking operator. The operator formed by a parallel architecture composed of an adaptive linear smoothing and an adaptive quadratic sharpening components, performs noise reduction and edge enhancement. The use of quadratic filters leads to a stronger enhancement of bright pixels than dark ones. Although this treatment is relevant for human perception, it could be detrimental for a post-processing. For instance, in seismic images black and white pixels represent positive and negative reflections, and it is important to enhance them similarly. In this case we propose to simply use a linear adaptive UM method. The proposed approach is developed in two steps. We first define the filter mask as a set of coefficients which quantifies the contribution of each pixel of the support in order to compute the output value. Depending on the type of images we want to process (natural, periodic or pseudo-periodic texture, seismic, etc.), we propose two different ways to compute these coefficients: locally adapted and recursively adapted. The second step consists in combining the pixel values according to this set of coefficients. The aim of this paper is to provide a robust enhancement filter in comparison with other 2D methods, and a 3D extension of this filter in order to process seismic images. The paper is organized as follows. In Section 2 the derivation of the best filter mask according to a performance criterion is presented. In Section 3, the contrast enhancement filter as a combination of a lowpass and a highpass filter is described. The two lowpass and highpass components are then developed and their properties discussed. In Section 4, results of image enhancement are presented in comparison to other approaches. Section 5 is devoted to the 3D extension of the filters, and its
S. Guillon et al. / Signal Processing 67 (1998) 237—254
239
application to the restoration of seismic images. Concluding remarks can be obtained in Section 6.
according to the minimization of a criterion defined later.
2. Filter mask
2.1. Local adaptive mask
Let us define the support S as the set of pixels likely to contribute in the filtering of the current pixel (x,y). Without loss of generality, we consider that S is a rectangular window of size K ]K : 1 2 K !1 K !1 S" (i, j)3 ! 1 ; 1 2 2
The simplest way to define the similarity between two pixels is to compare their grey levels. The approach is based on the observation that variations of gray levels inside a nontextured region are smaller than those between two different regions. Then, for a given location (i, j) of S, we estimate a continuity value D f (x#i, y#j)!f (x, y)D, and the coefficients m are computed by i,j m "C(D f (x#i,y#j)!f (x,y)D), (3) i,j where C is a positive decreasing function, such that C(d)P0 as d increases. A large number of such functions have been proposed in the bibliography for edge-preserving smoothing techniques [2]: inverse gradient filter [20], improved inverse gradient filter [19] and adaptive Gaussian weighted filter [16]. Hence a possible choice for C is
G C
C
D
K !1 K !1 ] ! 2 ; 2 2 2
DH
.
(1)
When scanning the image, the pixels included in the mask may belong to different regions. The principle of adaptive filtering then consists in adjusting the contribution of each pixel of S according to their similarity with the current pixel. For this purpose, we define the filter mask M as a set of coefficients m (in the range [0,1]) associated with i,j the support S: M"Mm 3[0,1]D(i, j)3SN. (2) i,j Each coefficient m is viewed as a level of confidence i,j for the pixel (x#i, y#j) to belong to the same region as the current pixel (x,y), i.e. m P1 for i,j pixels similar to the current pixel, respectively m P0 for the others. i,j In order to compute the set of coefficients M, a criterion must be defined to estimate the similarity between two pixels. A first criterion could be defined in comparing the two corresponding grey levels. The deviation between the two pixel values is calculated, and a value between 0 and 1 is assigned to m such that i,j m "1 when the deviation is zero, and m "0 for i,j i,j large deviation. This method is referred to a local adaptive approach. This estimation is fairly appropriate for natural images, for which each pixel to be filtered needs a specific mask. A second approach, devoted to periodic or pseudo-periodic textured images, is also presented. For this kind of images, the filter mask is stable or slow spatially-varying. Then we propose a recursive approach, in which the mask M evolves iteratively
A
B
( f (x#i,y#j)!f (x,y))2 m "exp ! , i,j 2d2
(4)
where d is a parameter controlling the width of the curve or the effective degree of similarity between the two positions. Its choice will be discussed in Section 4. The function has proved, through a large number of simulations, to give better results than those presented in [19,20]. Intuitively, from Eq. (4) we see that pixels belonging to the same region as the pixel (x,y) will have larger weighting coefficients than those located outside the region. Remark. In order to reduce the computation cost due to the Gaussian term in the above equations, note that the possible continuity values are limited to 255 values, for an 8-bit coded image. In practice a table can be used to estimate the m values. i,j 2.2. Recursive adaptive mask The idea of adapting recursively the filter mask is due to Salembier who derived an adaptive approach
240
S. Guillon et al. / Signal Processing 67 (1998) 237–254
for rank order-based filters [15] to process periodic or quasi-periodic textured images. We propose here another adaptation process of this method. The method, proposed by Salembier to define an adapted unweighted filter mask, consists in defining a search area, and assigning a weighting coefficient to each possible location. The current filter mask support is obtained by thresholding the set of coefficients — if a coefficient is greater than the threshold, the corresponding location is considered as belonging to the unweighted filter mask. An adaptation process is used to optimize (according to a criterion) the filter mask support by modifying the set of continuous values of the search area. In the approach we propose, a weighted filter mask is used instead of a binary one. All the pixels (x#i, y#j) belonging to a predefined support S, are taken into account according to their weight m — no arbitrary thresholding is necessary to i,j obtain a binary mask. The algorithm then operates in two steps: it first computes the coefficients m , i,j and secondly implements the linear filtering. The structure of the filter is given in Fig. 2 where the coefficients are adaptively modified by minimizing a prediction error. The predicted value p(x,y) we propose is a weighted mean over S*"SCM(0,0)N: f (x#i,y#j) + *m p(x,y)" (i,j)|S i,j . (5) m + (i,j)|S* i,j The optimality criterion we consider in this paper is the mean square error (MSE): J"E[(p(x,y)!f (x,y))2],
(6)
where E[)] is the mathematical expectation. Similar results are obtained with a mean absolute error (MAE) criterion. The minimization of the prediction error is performed with a steepest descent algorithm [21], chosen for the simplicity of the LMS algorithm, and its ability to deal with nonstationary signals. The minimization process is then reached by computing an estimate of the gradient of the criterion J with respect to the set of parameters M. If (mk ) i,j i,j|S is the current set of coefficients at step k, the next estimate at step k#1 is computed with the
Fig. 2. Filter structure.
following formula:
A
B
LJ mk`1"W mk !j , (7) i,j i,j Lm i,j where j is a parameter which controls the convergence of the algorithm, and W is a function maintaining the m values in the range [0,1] in order to i,j satisfy the definition of the filter mask M (Eq. (2)). The W function can be a simple clamping or an affine correction as follows: k !min k , mk`1" i,j i,j max !min k k where LJ k "mk !j , i,j i,j Lm i,j
(8)
(9)
min " min Mk N, k i,j (i,j)|S
(10)
max "max Mk N, (11) k i,j (i,j)|S In the derivation of the LMS algorithm, the gradient estimate of the criterion involved in Eq. (6) is computed by simply approximating the expectation with the instantaneous value. This leads to the following estimate: LJ LE[(p(x,y)!f (x,y))2] " Lm Lm i,j i,j Lp(x,y) "2E (p(x,y)!f (x,y)) Lm i,j Lp(x,y) K2(p(x,y)!f (x,y)) , Lm i,j
C
D (12)
S. Guillon et al. / Signal Processing 67 (1998) 237—254
241
with m !+ m f (x#r, y#s) Lp(x,y) f (x#i, y#j)+ p(x,y)!f (x#i, y#j) (r,s)|S* r,s (r,s)|S* r,s " "! . m )2 m (+ Lm + (r,s)|S* r,s i,j (r,s)|S* r,s So we obtain the following updating algorithm:
C
mk`1"t mk #2j(p(x,y)!f (x,y)) i,j i,j
A
BD
p(x,y)!f (x#i, y#j) ] . (14) mk + (r,s)|S* r,s As Eq. (14) is a nonlinear function of the coefficients m , we cannot easily carry out its convergence i,j properties. However, two well-known properties of the LMS algorithm [21] are still valid in this case: f the coefficients m exhibit an exponential coni,j vergence, with a time constant q inversely proportional to the parameter j: 1 qJ ; j f
(15)
the final misadjustment of the coefficients increases with j, and the mean square estimation error e is proportional to j: mi,j (16) E[e2 ]Jj. mi,j So the choice of j results from a compromise between the convergence speed (necessary in the case of non-stationary images) and the misadjustment of the coefficients m . i,j One can observe in the updating equation (14) that as p(x,y) is a predicted value of f (x,y), the set M evolves in order to improve the estimation p(x,y) of f (x,y). Indeed f if p(x,y) underestimates f (x,y) (i.e. p(x,y)( f (x,y)), then the values m for which f (x#i,y#j)' i,j p(x,y) are increased; f if p(x,y) overestimates f (x,y) (i.e. p(x,y)' f (x,y)), then the values m for which f (x#i,y#j)( i,j p(x,y) are increased. By copying this procedure, one can derive several updating equations which increase coefficients of pixels similar to the current pixel, and decrease other coefficients. For example, by introducing a similarity index r (in the range [!1,1]) between i,j the two pixels (x,y) and (x#i,y#j), which is max-
(13)
imum for two similar values and minimum for two different values, we propose the general following updating equation: (17) mk`1"W[mk #jr ]. i,j i,j i,j In the last section we will provide an example of such an updating function for the restoration of seismic images.
3. Contrast enhancement filter In order to improve the performances of the UM method, we propose to combine a lowpass and a highpass filter (see Fig. 3). The lowpass component g is used to smooth data in homogeneous areas, L while the highpass g is chosen as an edge detector H in order to sharpen details in the image. The filter output g(x,y) is given by g(x,y)"g (x,y)#ag (x,y), L H where
(18)
g (x,y)" + h f (x#i,y#j), L i,j (i,j)|S
(19)
g (x,y)" + w f (x#i,y#j). (20) H i,j (i,j)|S The coefficient a is a factor driving the edge enhancement effect. The operator we propose is then defined by the choice of the coefficients h and i,j w which are different combinations of the m i,j i,j
Fig. 3. Enhancement filter structure.
242
S. Guillon et al. / Signal Processing 67 (1998) 237–254
presented in the previous section. As reported in [10], to preserve the expected output of a uniform luminance input, the following conditions must hold:
S of pixels the most similar to the current pixel. 1 This threshold could be set arbitrary to a constant value (for example ¹"0.5), or be estimated adaptively by taking
+ h "1 and + w "0. (21) i,j i,j (i,j)|S (i,j)|S In the following we first give the expression of the lowpass and highpass filters, and finally propose an adaptive estimation of the scaling factor a.
1 ¹"mN " + m . (26) i,j Card(S) (i,j)|S With this value the class S is the set of pixels 1 having the highest levels of confidence.
3.1. Lowpass filter
3.2. Highpass filter: ¸aplacian like enhancement (¸¸E)
As the purpose of lowpass filtering is to reduce noise in the image without smoothing details or edges, we consider the following lowpass filter g (x,y): L g (x,y)" + h f (x#i, y#j ), L i,j (i,j)|S m i,j with h " . (22) i,j + m (k,l)|S k,l The condition (21) can be easily checked: m i,j "1. (23) + h " + i,j + m (i,j)|S (i,j)|S (k,l)|S k,l It results from the definition of the mask M, that g (x,y) is a weighted mean over the set of pixels the L most similar to (x,y) (pixels having high m values). i,j So the lowpass filter effectively performs noise reduction without edge smoothing. This filter belongs to the class of edge preserving smoothing filters (EPSF) widely studied in [2,14,16,19,20]. The performances of the lowpass filter can be improved by defining the following submask S : 1 S "M(i, j)3SDm '¹N, (24) 1 i,j where ¹3[0,1] is a threshold, and by performing the new lowpass filter: g1(x, y)" + h1 f (x#i,y#j), L i,j (i,j)|S1 m i,j with h1 " . (25) i,j + m (k,l)|S1 k,l According to the choice of the threshold ¹, g1(x,y) L becomes a weighted mean restricted over a class
In this section, we propose a new highpass filter defined by a linear combination of the coefficients of the filter mask M. We show that it is equivalent to a Laplacian filter, but that it exhibits an improved noise robustness. The highpass filter g (x,y) we propose is defined H by g (x,y)" + w f (x#i,y#j), H i,j (i,j)|S w "m !mN , i,j i,j with mN " 1 + m . C!3$(S) i,j (i,j)|S The set of coefficients W"(w ) verifies Eq. i,j + w " + (m !mN )"0. i,j i,j (i,j)|S (i,j)|S To point out the behaviour of the highpass ponent, we write g (x,y) as follows: H
G
(27) (21): (28) com-
g (x,y)" + (m !mN ) f (x#i,y#j) H i,j (i,j)|S " + (m !mN ) f (x#i,y#j) i,j (i,j)|S1 # + (m !mN ) f (x#i,y#j) i,j (i,j)|S2 " + Dm !mN D f (x#i,y#j) i,j (i,j)|S1 ! + Dm !mN D f (x#i,y#j), i,j (i,j)|S2
(29)
S. Guillon et al. / Signal Processing 67 (1998) 237—254
A
n(N!n) + f (x#i, y#j) (i,j)|S1 K N n
with S "M(i, j)3SDm *mN N 1 i,j set of pixels similar to current pixel,
B
f (x#i, y#j) + ! (i,j)|S2 N!n
S "M(i, j)3SDm (mN N 2 i,j set of pixels different from current pixel. So g (x,y) can be interpreted as the local difference H of pixel values belonging to two different regions — the ones belonging to S which are similar to the 1 current pixel, and the others belonging to S which 2 are different from the current pixel. The filter response can be studied in the following cases: f homogeneous area If the current pixel (x,y) is located on a quasi uniform area, then we can state that f (x#i,y#j)+f (x,y)Nm K1 ∀(i, j)3S i,j NmN K1
f
NDm !mN DK0 ∀(i, j)3S i,j Ng (x,y)+0. (30) H The output of the highpass filter is then negligible. presence of an edge If the pixel (x,y) is located across an edge (Fig. 4), the two classes S and S are characterized by 1 2 (i, j)3S Nm K1, 1 i,j (31) (i, j)3S Nm K0. 2 i,j So considering that Card(S)"N, Card(S )"n 1 and Card(S )"N!n, we have mN Kn/N and 2 n g (x,y)K + 1! f (x#i,y#j) H N (i,j)|S1 n ! + 0! f (x#i,y#j) N (i,j)|S2
A
A
B
B
243
n(N!n) K ( fM !fM ), @S1 @S2 N
(32)
with fM and fM the mean values of f estimated @S2 @S1 over S and S , respectively. Thus we easily show 1 2 that if (x,y) is located at the bottom of the edge, then g (x,y)(0 as fM (fM . If (x,y) is located at H @S1 @S2 the top of the edge, then g (x,y)'0 as fM 'fM . H @S1 @S2 And finally if (x,y) is located along the edge then g (x,y)"0 as fM "fM . H @S1 @S2 Globally the highpass filter is then equivalent to a classical Laplacian operator, but with much less noise sensitivity as we will see further (Section 4). More generally, g (x,y) is equivalent to a direcH tional second order derivative considered orthogonally to the edge. For example, a support S of size 3]3 located across an edge gives the following filter mask:
(33) with a a parameter which depends on the way the filter mask is implemented.
Fig. 4. Laplacian behavior.
244
S. Guillon et al. / Signal Processing 67 (1998) 237–254
Fig. 5. Filter structure with an adaptive a. Fig. 6. Adaptive a.
3.3. Scaling factor a As introduced in Eq. (18), a is a factor driving the enhancement effect. Choosing a constant value for a has the drawback of enhancing sharp edges more than smoothed ones because the enhancement is proportional to the sharpness of the detected edge. This implies an overshooting effect on the image (false black or white lines along edges), which could result in an unpleasant enhancement effect. This can be solved by applying a local adaptation of a depending on the sharpness of the edge. As the output of the highpass filter measures this sharpness, we use the value g (x,y) to control H a (Fig. 5): a"/(Dg (x,y)D). H
(34)
The enhancement function /()) must be continuous (/ : R`PR`) and satisfy f lim /(x)"0: a small g value is considered x?0 H as a false detection and no enhancement is required. f lim /(x)"0: a high g value corresponds x?= H to a sharp edge, and no enhancement is required. Lagrangian functions (Eq. (35)) satisfy all these constraints: a"/(Dg (x,y)D)"CDg (x,y)Dexp(!cDg (x,y)D). H H H
(35)
We can also use a simple piecewise linear approximation of this function (Fig. 6). Threshold ¹ controls noise removal and thre1 shold ¹ controls the overshooting effect. The 2 simulations presented in the last section are carried out with ¹ "5 and ¹ "50 (for 8-bit coded 1 2 images).
3.4. Iterative application of the enhancement filter The new proposed filter performs simultaneously edge enhancement and region smoothing. However, according to the level of noise and blur in the image, a single processing of the image could be not sufficient to achieve enough enhancement. We suggest in this case to apply iteratively the enhancement filter. Let g(x,y,n) be the filtered image at iteration n, then the filtering process becomes g(x,y,0)"f (x,y),
(36) g(x,y,n#1)"g (x,y,n)#ag (x,y,n) ∀n*0. L H This process is equivalent to that proposed in [14], but needs a smaller number of iterations — less than five iterations are sufficient in practice. We provide an example in the next section.
4. Application to image enhancement and filtering In this section we will describe some practical examples which demonstrate the performances of these new filters to restore blurred or noisy images. Simulations are performed with different classes of images — synthetic, natural and textured.
4.1. Application of the local filter In this section, the coefficients are calculated from the deviation between the current pixel value and the values of the pixels included in support S (Eq. (4)). The parameter d has to be set at first.
S. Guillon et al. / Signal Processing 67 (1998) 237—254
4.1.1. Choice of d In order to state how to choose the parameter d involved in the filter mask, we first study the behavior of the new filter on a noisy region with or without the presence of an edge. First, if the pixel (x,y) is located on a uniform area (mean value f ) corrupted by a white Gaussian 0 noise v(x,y) of variance p2: f (x,y)"f #v(x,y). (37) 0 Then the set of coefficients of the filter mask is defined by ( f (x,y)!f (x#i,y#j))2 m "exp! i,j 2d2 (v(x,y)!v(x#i,y#j))2 "exp! . 2d2
(38)
Assuming that v(x,y) is a Gaussian random variable, m is a realisation of a random variable m whose i,j density function p (m) can be expressed by m D m(D2~2)@2 d p (m)" m3[0,1] with D" . m p J2n J!ln m (39) Then, the mean mN and the variance p2 of these m coefficients can be obtained:
P
mN "E[m]" and
1
0
D mp (m)dm" m J2#D2
(40)
P
1 m2p (m)dm!mN 2 m 0 D D2 " ! . J4#D2 2#D2
p2 "E[(m!mN )2]" m
(41)
We can state from relations (40) and (41) that if we take D'2, then mN '0.8 and p2 )0.04, which m implies Nm KmN i,j
∀(i, j)3S
+ m ) g(x#i,y#j) Ng (x,y)" (i,j)|S i,j PB + m (i,j)|S i,j 1 K + g(x#i,y#j) Card(S) (i,j)|S
(42)
Ng (x,y)" + (m !mN ) ) g(x#i,y#j)K0. PH i,j (i,j)|S
245
This means that if d verifies d'2p, then the enhancement filter is equivalent to a smoothing filter, and the output variance is given by p2 p2" . g Card(S)
(43)
In the second case in which the pixel (x,y) is located across an edge, the choice of d depends on the height h of this edge. If a and b are the grey levels on the two sides of the edge, and if f (x,y)"a we can state that (i, j)3S Q f (x#i,y#j)Ka N m "m K1, 1 i,j 1 (i, j)3S Q f (x#i,y#j)Kb 2 (b!a)2 N m "m Kexp! i,j 2 2d2 h2 "exp! . (44) 2d2 The enhancement effect of the filter is effective only if the two classes S and S are well separated, i.e. if 1 2 m @m , which implies 2 1 (b!a)2 exp! @1 Q d@h. (45) 2d2 Finally, we can state that the choice of the parameter d must satisfy the following condition: 2p(d(h , (46) .*/ where p2 is the noise variance, and h the smallest .*/ height of edge we want to enhance. 4.1.2. Application In order to illustrate the behavior and the noise robustness of this enhancement filter, we compare the response of the UM and LLE methods. We first consider a blurred edge (h"110), and we compare results obtained by the highpass component: LLE (3]3, d"20) and Laplacian (3]3). In Fig. 7, one can observe that these two components are similar, and finally the edge enhancement is equivalent (Fig. 8). We then degrade the synthetic edge by adding a white Gaussian noise of variance p2"16. Results of the highpass components are provided in Fig. 9 — the support S is of size 5]5, and d"40 for the LLE. These results show the noise robustness of the LLE operator in comparison with the Laplacian.
246
S. Guillon et al. / Signal Processing 67 (1998) 237–254
Fig. 7. Comparison of the Laplacian behavior.
Fig. 8. Comparison of edge enhancement.
S. Guillon et al. / Signal Processing 67 (1998) 237—254
247
Fig. 9. Comparison of the Laplacian behavior on a noisy edge.
As expected, the new enhancement filter allows noise smoothing in uniform areas in the same time than edge enhancement (Fig. 10). In the second application, we propose to enhance the image of boat presented in Fig. 11. If we exercise the quadratic operators proposed by Mitra and Ramponi [8,10] on this image, we observe noise amplification in background area (such as the sky). The cubic operator proposed by Ramponi [9] is more robust against noise, but one can observe that some continuous edges are discarded. In Figs. 12 and 13 we provide the results obtained with the Mitra quadratic operator, the Ramponi cubic filter and the LLE operator, respectively. In these images, the contrast enhancement effect of the LLE operator is easily perceived without exhibiting noise amplification in background areas. An adaptive scaling factor a allows enhancement without overshoot — although the image obtained with an adaptive a is perceptually not as good as that obtained with a constant a, the absence of overshoot allows to use this image for a postprocessing.
4.2. Application of the recursive filter Some examples of enhancement of periodic or quasi-periodic textured images are presented to point out the interest of the recursive filter mask. Indeed, on such images, the optimal filter mask will be constant or slow spatially-varying. An adaptive evaluation of the filter mask M will then be more robust than a local one — if the pixel to be filtered is corrupted by noise, the coefficients of the filter mask will not be affected due to the high memory effect of M. The image to be filtered is the ‘canvas’ image presented in Fig. 14. If we process such an image by a quadratic or a cubic operator (Fig. 15: Ramponi’s cubic filter [9] with d"0.0003, Mitra’s quadratic filter [8] with a"512), we observe a noise amplification which gives unpleasant filtered images. On the contrary, results obtained with a 3]3 LLE filter appear more convincing in terms of contrast enhancement (Fig. 16: LLE method) — the texture is well enhanced without noise amplification. The image is scanned along columns so that the coefficients quickly converge.
248
S. Guillon et al. / Signal Processing 67 (1998) 237–254
Fig. 10. Comparison of edge enhancement on a noisy edge.
blurring. As the LLE operator combines a lowpass filter to remove noise and a highpass filter to restore edges, restoration with this operator will be more robust. We have iteratively filtered the noisy image with a"0.4, j"0.01 and S"[5]5] (Fig. 17). After four iterations the image becomes nearly binary with a good restoration of its main features, whose perception are difficult in the original image.
5. 3D extension and application to seismic images
Fig. 11. ‘Boat’ image.
In the second application, we propose to restore the ‘canvas’ image corrupted by a zero-mean Gaussian white noise with standard deviation p2"60. Classically additive Gaussian noise is removed by applying low pass filters, but at the expense of edge
The main purpose in seismic data analysis is the detection of ‘geological horizons’ separating homogeneous layers of rocks, sediments, etc. A 2D seismic image (Fig. 19) is composed of seismic traces, sinusoid-like waveforms, which are a record of reflected waves arising from impedance contrasts between strata. If a cycle is correlated laterally across many seismic traces on an image, it is called a seismic horizon. Filtering seismic images then consists in restoring these seismic horizons. As a seismic image is often
S. Guillon et al. / Signal Processing 67 (1998) 237—254
249
Fig. 12. Image enhancement.
Fig. 13. Local zoom.
K ]K ]K (corresponding support S). 1 2 3
G
C
D
K !1 K !1 S" (i,j,k)3 ! 1 ; 1 2 2
C C
D DH
K !1 K !1 ] ! 2 ; 2 2 2 K !1 K !1 ] ! 3 ; 3 2 2 M"Mm
,
3[0,1]D(i, j, k)3SN. i,j,k
(47) (48)
The output g(x,y,z) of the filter then becomes Fig. 14. ‘Canvas’ image.
extracted from a 3D bloc of data, the filtering process can be implemented over a 2D or 3D support. In the case of 3D processing, the filter mask M is naturally extended as a window of size
g(x,y,z)"g (x,y,z)#ag (x,y,z), L H
(49)
with the component g and g defined by a simple L H 3D extension: m f (x#i,y#j,z#k) i,j,k g (x,y,z)" + , L + m (i,j,k)|S i,j,k (i,j,k)|S
250
S. Guillon et al. / Signal Processing 67 (1998) 237–254
Fig. 15. Quadratic enhancement.
g (x,y,z)" + (m !mN ) f (x#i,y#j,z#k), H i,j,k (i,j,k)|S 1 with mN " + m . (50) i,j,k Card(S) (i,j,k)|S In the following, we propose local adaptive and recursively adaptive methods, specifically adapted to seismic images, due to the way of computing the set of coefficients M. To simplify the presentation, the equations are given in the 2D case, the 3D case being a simple extension by adding the third coordinate z.
5.1. Seismic filter mask In Section 2, we defined various criteria of similarity between two pixels in order to compute the coefficients m . Here we propose to improve the i,j similarity criterion for seismic images. Indeed two pixels could be said to be similar if they belong to the same seismic horizon, which means that their corresponding seismic traces are highly correlated. A seismic trace is a column of the image, and its representation is a sinusoid-like waveform (Fig. 18: intensity evolution along a column of the image). So for each position (x,y) we define the vector seismic trace by (x,y)"[ f (x!k,y)] . 53 k|*~N@2§N@2+
(51)
Fig. 16. LLE enhancement.
The similarity of two pixels is then locally estimated by measuring the correlation factor o of their two corresponding traces, which gives in the 2D case: 1 1 m " # o( (x,y), (x#i,y#j)) i,j 2 2 53 53 1 1 T (x,y) ) (x#i,y#j) 53 53 " # . 2 2 E (x,y)E ) E (x#i,y#j)E 53 53
(52)
S. Guillon et al. / Signal Processing 67 (1998) 237—254
251
Fig. 17. Iterative enhancement.
A second approach for computing these coefficients can be derived from the recursive scheme exposed in Section 2.2. Indeed we assume that the generic updating equation expressed by mk`1"t[mk #jr ], (53) i,j i,j i,j where r is a measurement of similarity betwen two pixels respecting r "#1, if (x,y) and (x#i,y#j) are similar, i,j r "!1, if (x,y) and (x#i,y#j) (54) i,j are widely different.
Fig. 18. Seismic trace.
Thus we could assume that m 3[0,1], and that i,j m "1 only if (x,y)"b (x#i,y#j) (b'0). As i,j 53 53 it takes into account the vicinity of the pixels, the estimation of the coefficients m is more robust to i,j noise than the one proposed in Section 2.1.
As the correlation factor o (Eq. (52)) of two seismic traces satisfies these conditions, we propose to define the updating equation by
A
B
T (x,y) ) (x#i,y#j) 53 mk`1"t mk #j 53 . i,j i,j E (x,y)E ) E (x#i,y#j)E 53 53
(55)
Thus we obtain an adaptive estimation of the set of coefficients M which tends to increase the weights
252
S. Guillon et al. / Signal Processing 67 (1998) 237–254
of the positions the most correlated to the current position. The parameter j drives the memory of the filter mask.
5.2. Application The two proposed schemes we have presented are used to process the seismic image shown in Fig. 19. We first process this image with a 3D locally adapted mask, whose weights are evaluated with the correlation factor of the seismic traces (Eq. (52)), S"[5]5]5], a seismic trace of size 11 (N"5) and an enhancement factor a"0.1. The result is presented in Fig. 20. To evaluate the filtering effect of our approach on this class of images, we should focus our attention on the main features of these images: seismic hor-
izons, channels (located in the upper left of the image) and faults (right side of the image). Important features for expert interpretation are preserved while an efficient cancellation of noise is performed. We compare this result with the filtering using a recursive mask, where the coefficients evolve according to the correlation factor of the seismic traces (Eq. (55)). We scan the image in the column direction so that the filter mask is slow spatiallyvarying. The image is processed with an S"[5]5]5], a vector seismic trace of size 11, an updating coefficient j"0.01, and an enhancement factor a"0.3. The result is given in Fig. 21. The evaluation of the result is still done by focusing our attention to the features of this image: seismic horizons appear more regular than those resulting from the local filter (Fig. 20), but the channel and the fault are not as well preserved as with the local filter. Finally, the two approaches (local and recursive) can be combined if we are able to detect the position of the two main features (fault and channel). If we note d(x,y) the probability function of such a detection (d(x,y)P1 if a feature is detected), then the output image is obtained by
Fig. 19. Seismic image.
g(x,y)"g (x,y) ) d(x,y)#g (x,y) ) (1!d(x,y)), (56) -0# 3%# with g and g respectively being the local and -0# 3%# recursive filtering. From Fig. 22, one can observe that the main features are well preserved, and horizons regularized. This method combines the advantages of the two previous ones.
Fig. 20. Seismic filtered image with a local mask.
Fig. 21. Seismic filtered image with a recursive mask.
S. Guillon et al. / Signal Processing 67 (1998) 237—254
253
Notations f ()) g()) (x,y) (x,y,z) S (i,j) (i,j,k) Fig. 22. Seismic filtered image combining the local and recursive mask.
M"(m ) (or (m )) i,j i,j,k ()M ) E[)]
The detection d(x,y) we used is based on a measure of dispersion of local gradients [4] not described in this paper.
6. Conclusion In this paper new nonlinear image enhancement filters have been presented. They are based on the choice of a locally-adapted filter mask or a recursive adaptive process for computing the mask coefficients. The filter mask defines the set of surrounding pixels most closely correlated to the current pixel. Then two new enhancement techniques, called Laplacian like enhancement (LLE1 and LLE2), have been introduced as an improvement of the classical unsharp masking (UM) method. We have shown, through tests on synthetic and real images, the improvements brought by these algorithms: noise reduction and edge enhancement are performed simultaneously. Comparisons with other techniques have also been provided. Finally, these new approaches, including a 3D extension, have been applied to practical applications dealing with seismic images restoration. From a computational point of view, the use of a locally-adapted filter mask is not more costly than other UM algorithms (parallel implementation is possible), while the adaptive process for computing the coefficients leads to a more costly algorithm.
input image output image position of a pixel on an image position of a voxel on a 3D bloc of data filtering support position of a pixel on a 2D support position of a voxel on a 3D support set of levels of confidence (respectively 2D or 3D) mean operator mathematical expectation
Acknowledgements The authors wish to thank ELF Aquitaine for having provided a large number of seismic images, and for partially supporting this work. They are also indebted to the academic partners of the joint project CNRS — ELF Aquitaine. They also thank the anonymous reviewers for their helpful comments that aided a lot in making this paper more clear and complete.
References [1] R.J.P. de Figueiredo, S.C. Matz, Exponential nonlinear Volterra filters for contrast sharpening in noisy images, in: Proc. ICASSP’96, Vol. 4, pp. 2263—2266. [2] J.M.H. Du Buf, T.G. Campbell, A quantitative comparison of edge-preserving smoothing techniques, Signal Processing 21 (4) (1990) 289—301. [3] S. Guillon, P. Baylou, New noncausal adaptive quadratic filters for image filtering and contrast enhancement, in: Proc. EUSIPCO’96, Trieste, Italy, Vol. 1, 1996, pp. 583—586. [4] S. Guillon, P. Baylou, N. Keskes, Detection and/or determination method of features related to singular points, Submitted for a French patent No. 9708602, INPI (Institut National de la Protection Industrielle), July 1997. [5] S. Guillon, P. Baylou, M. Najim, Robust nonlinear contrast enhancement filters, in: Proc. ICIP’96, Lausanne, Switzerland, Vol. 1, 1996, pp. 757—760.
254
S. Guillon et al. / Signal Processing 67 (1998) 237–254
[6] A.K. Jain, Fundamentals of Digital Image Processing, Prentice-Hall, Englewood Cliffs, NJ, 1989. [7] J.S. Lee, Digital image enhancement and noise filtering by use of local statistics, IEEE Trans. Pattern Anal. Machine Intell. 2 (1980) 165—168. [8] S.K. Mitra, H. Li, I. Lin, T. Yu, A new class of nonlinear filters for image enhancement, Proc. ICASSP’91, pp. 2525—2528. [9] G. Ramponi, A simple cubic operator for sharpening an image, Proc. IEEE Workshop on NSIP, Thessaloniki, Greece, 20—22 June 1995, pp. 963—966. [10] G. Ramponi, P. Fontanot, Enhancing document images with a quadratic filter, Signal Processing 33 (1) (1993) 23—24. [11] G. Ramponi, G.L. Sicuranza, Quadratic digital filters for image processing, IEEE Trans. ASSP 36 (6) (June 1988) 1263—1285. [12] G. Ramponi, G.L. Sicuranza, Image sharpening using a polynomial filter, Proc. ECCTD’93, Davos, Switzerland, pp. 1431—1436. [13] G. Ramponi, N. Strobel, S. Mitra, T.H. Yu, Nonlinear unsharp masking methods for image contrast enhancement, Journal of Electronic Imaging 5 (3) (July 1996) 353—366.
[14] P. Saint-Marc, J.S. Chen, G. Medioni, Adaptive smoothing: a general tool for early vision, IEEE Trans. Pattern Anal. Machine Intell. 13 (6) (June 1991) 514—529. [15] P. Salembier, Adaptive rank order based filters, Signal Processing 27 (1) (1992) 1—25. [16] M. Spann, A. Nieminen, Adaptive gaussian weighted filtering for image segmentation, Pattern Recognition Letters 8 (1988) 251—255. [17] S. Thurnhofer, Quadratic Volterra filters for edge enhancement and their application in image processing, Ph.D. Thesis, University of California, Santa Barbara (1994). [18] A. Vanzo, G. Ramponi, G.L. Sicuranza, An image enhancement technique using polynomial filters, in: Proc. ICIP’94, Austin, TX, 13—16 November 1994. [19] D. Wang, O. Wang, A weighted averaging method for image smoothing, Proc. IEEE Conf. Com. Vision and Pattern Recognition, Miami Beach, Florida, USA, 22—26 June 1986, pp. 981—983. [20] D. Wang, A. Vagnucci, C.C. Li, Gradient inverse weighted smoothing scheme and the evaluation of its performances, CGIP 15 (1981) 167—181. [21] B. Widrow, S.D. Stearns, Adaptive Signal Processing, Prentice-Hall, Englewood Cliffs, NJ, 1985.