- Email: [email protected]
Image segmentation by contours and regions cooperation
Signal Processing 78 (1999) 329}347
Image segmentation by contours and regions cooperation A. Chehikian* Laboratoire des Images et des Signaux, Institut National Polytechnique de Grenoble and Universite& Joseph Fourier Grenoble, INPG, 46 Avenue Felix Viallet, 38031 Grenoble Cedex, France Received 13 July 1998; received in revised form 1 February 1999
Abstract In this paper we propose an algorithm for image segmentation in which contours and regions cooperate. In order for this cooperation to be e$cient, the role of contour pixels and region pixels should be considered di!erently. Region pixels carry information about the content of the image. They are de"ned on a lattice identical to the image lattice. Contour pixels carry information about regions boundary. They should be located between region pixels, therefore on a lattice distinct than the image lattice. A way to make these two statements concord is to achieve segmentation on a double-sized lattice where sites of both even coordinates are reserved for image pixels, whereas the remaining sites are reserved for contour pixels. We therefore propose a gradient operator, based on B-spline interpolation, which locates contour pixels between image pixels. Contours detected in this manner are used in an image segmentation algorithm which makes contours and region cooperate. 1999 Elsevier Science B.V. All rights reserved. Zusammenfassung In diesem Beitrag schlagen wir einen Algorithmus zur Bildsegmentierung vor, in welchem Konturen und Bereiche kooperieren. Um diese Kooperation e$zient zu gestalten, sollten die Konturpunkte und Bereichspunkte unterschiedlich behandelt werden. Bereichspunkte beinhalten Informationen uK ber den Bildinhalt. Sie sind de"niert uK ber eine Gittergleichheit zum Bildgitter. Konturpunkte beinhalten Informationen uK ber Bereichsgrenzen. Sie sollten zwischen Bereichspunkten liegen und damit auf einen vom Bildgitter unterschiedlichen Gitter. Ein Weg, diese Darstellungen in Einklang zu bringen, ist die Scha!ung der Segmentierung auf einem Gitter doppelter GroK {e, wobei PlaK tze beider geraden Koordinaten fuK r Bildpunkte und die uK brigen PlaK tze fuK r Konturpunkte reserviert sind. Wir schlagen hierzu einen Gradientenoperator basierend auf B-Spline-Interpolation vor, welcher die Konturpunkte zwischen den Bildpunkten lokalisiert. Auf diese Weise detektierte Kontouren werden in einem Bildsegmentierungsalgorithmus benutzt, welcher Konturen und Bereiche kooperieren laK sst. 1999 Elsevier Science B.V. All rights reserved. Re2 sume2 On propose dans cet article un algorithme de segmentation d'image par coopeH ration entre contours et reH gions. Pour que cette coopeH ration soit e$cace, les pixels de contour et ceux de reH gions doivent e( tre consideH reH s di!eH remment. Les pixels formant une reH gion apportent une information sur le contenu de l'image et sont donc deH "nis sur une grille d'eH chantillonnage identique a` celle de l'image. Les pixels deH "nissant les contours apportent, eux, une information sur la frontie`re des reH gions. LocaliseH s entre les reH gions, ils doivent e( tre deH "nis sur une grille d'eH chantillonnage distincte de celle
* Corresponding author. Tel.: #334-7657-4358; fax: #334-7657-4790. E-mail address: [email protected] (A. Chehikian) 0165-1684/99/$ - see front matter 1999 Elsevier Science B.V. All rights reserved. PII: S 0 1 6 5 - 1 6 8 4 ( 9 9 ) 0 0 0 7 3 - 0
330
A. Chehikian / Signal Processing 78 (1999) 329}347
de l'image. Pour assurer la coheH rence de chaque type de pixel, on propose de reH aliser la segmentation sur une grille d'eH chantillonnage de reH solution double de celle de l'image. Sur cette grille, les sites de coordonneH es toutes deux paires sont reH serveH s aux pixels d'image, les autres sont reH serveH s aux pixels de contour. Nous proposons un opeH rateur d'estimation du gradient, baseH sur l'interpolation B-spline, qui situe les contours entre les pixels d'image. Nous utilisons les contours ainsi deH tecteH s, dans un algorithme de segmentation mettant en wuvre une coopeH ration contours-reH gions. 1999 Elsevier Science B.V. All rights reserved. Keywords: Image segmentation; Contour detection; B-spline interpolation
1. Introduction Image segmentation is a necessary step in image understanding or low bit-rate coding. The goal of image segmentation is to split the image into regions. A region is a set of connected pixels having roughly the same value for some property, usually the gray level (or some other scalar property), surrounded by a closed contour representing its boundary. Two classes of image segmentation algorithms may be found in the literature. In the "rst, closed contours are detected, and regions are then deduced from their boundary by "lling the interior with the same label. Algorithms from the second class split the image into regions using some homogeneity criterion. Boundaries are then deduced if needed. Both type of algorithms have advantages and drawbacks. Contour detection is sensitive to noise and may result in both over-segmentation and under-segmentation. On the other hand, algorithms for contour detection are easy to implement because of the local nature of operations. Region based algorithms are usually more demanding in computation while being less sensitive to noise. Contour and region segmentation processes are usually considered as dual. One may take advantage of the two methods to obtain the best segmentation according to some segmentation quality criterion. However, such a cooperation requires that the role of each kind of pixel is well de"ned. Region pixels carry informations about an image property, they are therefore de"ned on the image lattice. Contour pixels carry informations about region boundaries, they should be located between regions and therefore should be de"ned on a distinct lattice.
In this paper, we propose an algorithm which makes contours and regions cooperate to obtain fast and e$cient image segmentation. In this algorithm, contours are detected "rst, regions are then obtained by "lling the interior of closed contours with the same label. The resulting regions are then used to validate signi"cant contours: close contours that separate two regions which do not satisfy some merging criterion. Final segmentation results from this cooperative process. With this approach, contour detection is a key process for which we propose an e$cient algorithm. Most contour detection algorithms are based on the evaluation of spatial gradient or Laplacian. Pixels are labeled contour pixels if an extrema of the gradient or a zero crossing of the Laplacian is observed. Such contour detection schemes raise the question of evaluating derivatives from a sampled signal, and also the question of the in#uence of noise when computing such derivatives. Canny [2], in his well-known approach, considers the case of a continuous edge modeled by a unit step. He looks for a "lter which exhibits a maximum output when such a step is present. To be optimal, such a "lter should: maximize the probability of detecting a contour, minimize the probability of detecting false contours due to noise, minimize the average displacement of the contour in the presence of noise and "nally, maximize the average distance between the detected false contours. The optimal "lter he proposes, is approximated by the "rst derivative of a Gaussian, so that contour detection reduces to searching the maximum of the "rst derivative of the output of a smoothing "lter. Deriche [4], considering an in"nite impulse response "lter, gives an exact solution for Canny's problem. Shen and Castan [3] also consider the case of a continuous step as
A. Chehikian / Signal Processing 78 (1999) 329}347
a model of an edge. They look for a "lter having the following properties: when only the step is present, the energy of the "rst derivative of the output should be maximum, while this energy should be minimum when only the noise is present. They conclude that the signal must be smoothed by a "rst-order in"nite impulse response "lter, the contour is detected when the "rst derivative of the output exhibits an extremum. It is interesting to note that all of these approaches consider the case of a continuous signal; the solutions that are proposed are extended to the discrete case by sampling the continuous impulse responses. In the contour detection algorithm we propose, we consider the more realistic case of a contour modeled by a discrete step. As we aim to estimate the "rst and/or the second derivative in the presence of such a signal, not only on sample points but between them, and because it is computationally e$cient, we have chosen to "t a continuous cubic spline to the discrete samples. However, in the presence of step contours, spline functions tend to oscillate and therefore provide multiple extrema around the true contour. It is therefore necessary to regularize the sampled signal in order to ensure the uniqueness of the gradient maximum or the Laplacian zero-crossing, while keeping a good localization of the contour. This paper is organized as follows. In Section 1, we remind some fundamental results concerning spline interpolation. These results are extracted from [5,6], readers are suggested to refer to these papers for more details. In Section 2, we consider the approximation of a discrete step by a cubic spline and deduce the values of the "rst and second derivatives at sample points. This will point out that cubic splines are not convenient for step contour detection, unless the bandwidth be reduced by an appropriate "lter. We then propose a set of regularization "lters which ensure the uniqueness of the gradient maximum and the zero-crossing of the Laplacian. We deduce the resulting 2D gradient operator, and show that a strong similarity exists between this operator and those of Shen}Castan and Deriche. However, if the continuous response underlying the discrete response of the proposed
331
operator is up-sampled on a double-sized lattice, the resulting discrete operator performs signi"cantly better than these popular edge detectors. In Section 3, we derive the interpolated operator which provides contour pixels located between the image pixels. We de"ne a detection criterion which makes more sense, in the discrete case, than the localization criterion, K, de"ned by Canny. We show that the best smoothing "lter, regarding the detection behavior, is again the one which induces the lowest computation cost. We illustrate on some examples, the bene"t we can get from our operator comparing with the well-known Deriche or Shen}Castan operators. In Section 4, "nally, we propose a very simple algorithm for image segmentation which makes contours and regions cooperate, and give some experimental results.
2. Discrete signal approximation by a continuous spline: some fundamental relations Unser et al. [5,6] consider the problem of approximating a signal discretized on a regularly spaced grid, by a spline function of degree n. Such a spline function is continuous and equal to a polynomial of degree n on each interval [k, k#1] if n is odd, on [k!, k#] if n is even (k3Z), continu ously derivable up to the order n!1, it is entirely de"ned by the set of coe$cients +y(k)3l , by the equation > gL(x)" y(k) bL(x!k), x3R, I\
(1)
where bL(x) is the continuous B-spline of order n. Such a function will exactly approximate a discrete signal +s(k), if +gL(k),"gL(x)" "+y(k),*+bL(k),"+s(k),, VI
(2)
where +bL(k), is the discrete B-spline. From Eq. (2), and considering the z-transform, we can write Eqs. (3) and (4) which de"ne a reversible transformation. Eq. (3) de"nes the direct transform giving the set of coe$cients +y(k), from the set
332
A. Chehikian / Signal Processing 78 (1999) 329}347
+s(k), whereas Eq. (4) de"nes the inverse transform giving +s(k), from the set +y(k),.
at a higher sampling rate through the use of the inverse transform:
X y(z)"s(z)BL(z)\, +y(k),&
(3)
X s(z)"y(z)BL(z). +s(k),&
(4)
+s (k),"+gL (k),"+[y(k)] ,*+bL (k),, (8) K K tK K where m is the zooming rate, +[y(k)] , is the set tK +y(k), expanded by the factor m:
In these equations, BL(z)\ is the transfer function of the direct B-spline "lter of order n, whereas BL(z) is the transfer function of the indirect B-spline "lter. A useful property of these transforms lies in the fact that the derivatives of the continuous spline being continuous splines (of lower degree) they are fully characterized by their coe$cients which can be easily deduced from the coe$cients y(k). First and second derivative values at sample points are obtained by
*gL(x) *x
and +bL (k), is the up-sampled discrete B-spline. K Eqs. (3), (5) and (6) will be used to obtain the value of the "rst and second derivative at sample points. In addition, Eq. (8) will provide the same value on a grid which is m times denser. To use these relations we need to know the coe$cients of the direct and indirect "lters. They are obtained from the recursive equations:
(10)
((n#2)/2#k/m)bL\(k#m)#(n/2!k/m)bL\(k) K K , cL (k)" K n
(11)
(5)
"+y(k#1)!2y(k) VI #y(k!1),*+bL\(k),,
with the starting conditions: m m 1 for ! )k) , 2 2 b (k)" K 0 else,
(6)
1 for 1!m)k)0, c (k)" K 0 else.
(7)
3. Edge detection using the cubic spline transform
where +cL\(k), is the discrete shifted B-spline: cL\(k)"bL\(k#).
(9)
((n#1)/2#k/m)cL\(k)#((n#1)/2!k/m)cL\(k!m) K K bL (k)" , K n
"+y(k)!y(k!1),*+cL\(k),, VI
*gL(x) *x
y(k) if k "mk, +[y(k)] ,"y(k )" tK 0 elsewhere,
These relations should be interpreted as follows. The coe$cients of the "rst derivative are obtained by the di!erence [y(k)!y(k!1)], values of the "rst derivative at sample points are further obtained from the convolution with the shifted Bspline of degree n!1. In the same way, the coe$cients of the second derivative are obtained by [y(k#1)!2y(k)#y(k!1)], samples of the second derivative are then obtained by convolving with the symmetric B-spline of degree n!2. Another useful property of the spline transform lies in the fact that the signal may be reconstructed
3.1. One-dimensional formulation In the following, we consider a discrete 1D step contour modeled by a discrete unit step +u(k),:
u(k)"
1 for k*0,
(12)
0 else.
We shall also consider the unit impulse +i(k),:
i(k)"u(k)!u(k!1)"
1 if k"0, 0 else.
(13)
A. Chehikian / Signal Processing 78 (1999) 329}347
As we are only interested in the "rst and second derivative at sample points, a cubic spline transform is quite appropriate and results in lower computation cost than if a higher degree spline transforms were used. In order to compute the values of the "rst and second derivatives, we need to know the transfer functions of the direct spline "lter of order 3 to obtain the spline coe$cients, the indirect shifted "lter of order 2 to obtain the sampled "rst derivative, and the indirect "lter of order 1 to obtain the sampled second derivative. Parameters of these transfer function are obtained from Eqs. (10) and (11): B(z)"1,
z#1 , C(z)" 2
One can notice that if indirect "lters are of "nite response, the direct spline "lter is of in"nite impulse response. Writing 6 B(z)\" z#4#z\
1 1!a 1 " # !1 , 1#a 1!az 1!az\ where a"(3!2(0, it is easy to show that the discrete impulse response +r(k), of the direct spline "lter of order 3 is given by 1!a r(k)" a I . 1#a
(14)
Eq. (3) yields +y(k),"+u(k),*+r(k),.
(15)
It results from Eqs. (5) and (6) that the spline coe$cients of the "rst and second derivatives are given by Eqs. (16) and (17), respectively: +y(k)!y(k!1),"+u(k)!u(k!1),*+r(k), "+i(k),*+r(k),"+r(k),,
+y(k#1)!2y(k)#y(k!1),"+r(k#1)!r(k),. (17) Therefore, the values of the "rst and second derivatives at sample points are given by
+u (k),"
r(k)#r(k#1) 2
(18)
(19)
1!a aI for k*0, 2 Nu (k)" 1!a a\I> for k(0, 2 +u (k),"+r(k#1)!r(k),
6 B(z)\" . z#4#z\
333
(16)
(1!a) ! aI for k*0, (1#a) Nu (k)" (1!a) a\I> for k(0. # (1#a)
Fig. 1(a) displays the continuous cubic spline approximating the discrete 1D step, as well as its "rst derivative, whereas Fig. 1(b) displays the second derivative. Sample points are marked with `a. We can see on these "gures that, in the vicinity of the step, the maximum of the "rst derivative and the zero crossing of the second derivative occur exactly at the midpoint of the gray-level transition, i.e. between the two image pixels. However, these derivatives exhibit multiple extrema and zero-crossings, respectively. Cubic B-spline interpolation cannot be used in this way for contour detection unless the signal bandwidth be reduced by an appropriate smoothing "lter.
3.2. Image bandwith reduction for contour detection by cubic spline transform We need to de"ne a smoothing "lter of impulse response +F (k),, which ensures that the second N derivative of the cubic spline which interpolates the discrete step smoothed by this "lter, will exhibit a unique zero crossing. Looking for a "lter of
334
A. Chehikian / Signal Processing 78 (1999) 329}347
Fig. 1. (a) First and (b) second derivative of the cubic spline approximating a step edge.
E The zero crossing of the Laplacian should be, for symmetry reasons, located at x"!0.5:
transfer function of the form X F (z)"B(z) f (z) +F (k),& N N N
u( (!1)"!u( (0)Nf (1)"f (!1). N N
will simplify the problem. Indeed, the coe$cients of the cubic spline interpolating the smoothed signal, will be given by X y( (z)"s(z)F (z)B(z)\"s(z) f (z) +y( (k),& N N
E The zero crossing of the Laplacian should be unique: u( (k))0, ∀k*0Nf (k#1)!f (k))0, ∀k*0, N N u( (k)*0, ∀k(0Nf (k#1)!f (k)*0, ∀k(0. N N (22)
X +s(k), +f (k), & * N and the values of the "rst and second derivative of the smoothed step +u( (k), at sample points will be, from Eqs. (18) and (19): f (k#1)#f (k) N , u( (k)" N 2
These three conditions are satis"ed by any "lter of positive, symmetric, unimodal impulse response: a Gaussian "lter for example, would be a good candidate. However, because of their interesting properties, we focus on the family of "lters of transfer function:
N 1!b 1!b f (z)" , N 1!bz 1!bz\
u( (k)"f (k#1)!f (k). N N
b"e\?, a'0Nb3]0, 1[ ,
We can now express some conditions that are to be satis"ed.
whose impulse response are of the form
E The gradient estimated on a positive step contour, should be positive:
+ f (k)," N
f (k)*0, N
∀k.
(21)
(20)
1!b N "k"N\ N\ b I # c "k"G G 1#b (p!1)! G
c '0 ∀i. G
(23)
,
(24)
A. Chehikian / Signal Processing 78 (1999) 329}347
Such impulse responses are positive, symmetric, and unimodal and therefore satisfy conditions (20)}(22). The "rst two members of this family of "lters, as it will be seen in the following, lead to a very popular contour detection operators. Their transfer function and impulse response are given by 1!b 1!b f (z)" , 1!bz 1!bz\
1!b b I , + f (k)," 1#b
(25)
(26)
335
1!2b#b 1!2b#b , f (z)" 1!2bz\#bz\ 1!2bz#bz
+ f (k),"
1#b 1!b b I "k"# 1#b 1!b
.
(27)
(28)
Fig. 2 displays the continuous cubic spline which approximates the discrete step smoothed by the "lter F , along with its "rst derivative (Fig. 2(a)) and its second derivative (Fig. 2(b)). Fig. 3 displays the same results when using the "lter F instead of F . In both "gures, a is set to 2. One can observe from these "gures that the two "lters, even of di!erent order, provide continuous
Fig. 2. (a) First and (b) second derivative of the cubic spline approximating a step edge smoothed by F .
Fig. 3. (a) First and (b) second derivative of the cubic spline approximating a step edge smoothed by F .
336
A. Chehikian / Signal Processing 78 (1999) 329}347
splines and spline derivatives of the same shape: piecewise cubic for the spline, piecewise quadratic for its "rst derivative and piecewise linear for its second derivative. One may wonder if the use of a high-order smoothing "lter, which of course induces a higher computation cost, is more bene"cial than the use of a lower-order smoothing "lter. This point will be discussed later. We would like to "rst de"ne the 2D operators resulting from the use of these smoothing "lters. 3.3. Deriving a 2D gradient operator When using the smoothing "lters de"ned in the previous section, it results from relations (3), (5) and (4) that the gradient values at sample points are obtained by
whereas the transfer functions of the Shen}Castan "lters may be expressed by
1!b 1!b , F (z)" 1 1!bz 1!bz\
One can see that the derivation "lters di!er only by a normalizing factor, whereas the smoothing "lters di!er by the indirect spline "lter which is not present in the Shen}Castan's operator. Using the smoothing "lter of order 2: z#4#z\ F (z)" 6 ;
(36)
;
where (z!z\) H (z)"f (z) N N 2
s( (z , z )"s(z , z )H (z )F (z ), (30) N N s( (z , z )"s(z , z )H (z )F (z ). (31) N N These relations are very similar to those proposed by Shen}Castan and Deriche for the estimation of the gradient vector components. It is worthwhile to compare the proposed operator to theirs. Using the smoothing "lter of order 1:
z#4#z\ 1!b 1!b F (z)" , 6 1!bz 1!bz\
1!2b#b 1!2b#b , 1!2bz#bz 1!2bz\#bz\ (37)
(29)
and p is the order of the smoothing "lter. In the 2D case, smoothing is performed in the two orthogonal directions, whereas derivation is performed only in the direction of the gradient's component. It results that the two gradient components will be obtained by
1!2b#b 1!2b#b , 1!2bz#bz 1!2bz\#bz\
(z!z\) H (z)" 2
(z!z\) s( (z)"s(z) f (z) "s(z) H (z), N N 2
(z!z\) 1!b 1!b H (z)" , 2 1!bz 1!bz\
b(1#b) 1!b 1!b H (z)" (z!z\) . (35) 1 (1!b) 1!bz 1!bz\
X s( (z)"[[s(z) F (z)]B(z)\] (1!z\) (1#z), s( (k)& N 2
(34)
whereas the transfer functions of the Deriche "lters (case u equal to 0) may be expressed by F (z)"(c z#c #c z\) " 1!2b#b 1!2b#b ; , 1!2bz#bz 1!2bz\#bz\
(38) 1#b H (z)" (z!z\) " 1!b ;
1!2b#b 1!2b#b , 1!2bz#bz 1!2bz\#bz\ (39)
where
(32)
b (b\!4a!b) c " , (1!b) (b\#2a!b)
(33)
b ((a!1)b\#(a#1)b) c " . (1!b) (b\#2a!b)
A. Chehikian / Signal Processing 78 (1999) 329}347
Once more, the derivation "lters di!er only by a normalizing factor, whereas the smoothing "lters are very similar. Moreover, it can be shown that when a goes to zero, the coe$cients c and c in the Deriche smoothing "lter go to 4/6 and 1/6, respectively (a equal to 1 yields c and c equal to 0.6886 and 0.1557, respectively). One may conclude that our proposed operator and Deriche's operator, although developed from two very di!erent approaches, are identical when smoothing becomes e!ective. However, we must repeat that these comparisons make sense only when considering the discrete "lters, as the continuous "lters are quite di!erent. We compare in the next section, the response of Deriche and Shen}Castan "lters to those of the operators we propose.
3.4. Comparing with Deriche and Shen}Castan xlters A common characteristic of these "lters is that they are de"ned to optimize some criterions when considering a contour modeled by a continuous step. The results which are obtained are then extended to the discrete case by sampling the continuous impulse response. Figs. 4 and 5 show the response of the continuous and discrete "lters when applied to a step contour without and with an additive noise. The smoothing parameter of each "lter has been chosen to ensure the same signal-tonoise criterion R"1, as de"ned in Canny's approach. In both "gures, the ratio A/n of the step amplitude to the noise standard deviation is 8. One may observe in Fig. 4 illustrating the continuous case that, for both "lters, even in the presence of noise, the maximum of the response is quite well localized, due to the moderate noise environment. Considering the discrete case, illustrated in Fig. 5, results in two consequences. First, due to the violation of the sampling theorem, the exact position of the contour is lost. Indeed, any continuous step localized at x 3]!1, 0] will be described by the same discrete step. Second, the response of the discrete "lters exhibits a maximum in form of a plateau, when applied to an un-noisy step, whereas it exhibits a true maximum, but randomly
337
located at x"0 or !1, when applied to a moderately noisy step. As a matter of fact, extending Canny's localization criterion, initially de"ned from a continuous model of contour, to the discrete case is at least improper. Fig. 6 shows the normalized response of our proposed "rst-order "lter. Because in our approach the contour we consider is modeled by a discrete step, we only consider the response of the "lter to such a discrete step. The sets of points marked with `a constitute the response of the "lter discretized on the input lattice. When considering these responses, one may conclude, in accordance with Eqs. (33) and (35), that the proposed "lter exhibits the same behavior as the one of Shen}Castan (the same conclusion would be derived from a comparison between our second-order "lter and the one of Deriche). However, because we interpolate the discrete step by a continuous spline, we are able to provide the underlying continuous response and the up-sampled discrete response as well. The sets of points marked with ` 6 *a constitute the discrete response up-sampled on a double sized lattice. Such response exhibit a true maximum at x"!0.5, the position of which is very less a!ected by the noise than the response of the "lter discretized on the input lattice. The following section de"nes an interpolated operator for contour detection which reduces the uncertainty on contour position and provides contour sites located between the image pixels, thus facilitating image segmentation.
4. An interpolated operator for contour detection The structure of such an operator is very similar to the one of the previous 2D operator except that it includes an interpolating step according to Eqs. (8) and (9). For the sake of simplicity, we de"ne "rst the operations to perform to interpolate without derivation and next when interpolating with the "rst-order derivation. From Eq. (8), we can write +s( (k),"+[y( (k)] ,*+b(k), t t X ' & s( (z )"y( (z ) B(z ), t
(40)
338
A. Chehikian / Signal Processing 78 (1999) 329}347
Fig. 4. Shen}Castan and Deriche "lters response to a continuous step.
A. Chehikian / Signal Processing 78 (1999) 329}347
Fig. 5. Shen}Castan and Deriche "lters response to a discrete step.
339
340
A. Chehikian / Signal Processing 78 (1999) 329}347
Fig. 6. The proposed "rst order "lter's response to a discrete step.
with from Eq. (10), 0.125z \#z\#2.875z \#4#2.875z #z #0.125z B(z )" 6 which can be factorized into z \#4#z z\#4z \#6#4z #z B(z )" 6 8
(1!z \) C(z )"(1!z \) can be factorized into
(41)
(1!z \)C(z )
where one may recognize (up to a multiplicative factor 2) in N(z ) the binomial "lter of order 4. In the same manner, from Eqs. (5) and (8):
"(z !z \)
"B(z ) N(z ),
+s( (k),"+[y( (k)!y( (k!1)] ,*+c(k),, t t s( (z )"y( (z ) (1!z \) c(z ), t where, from Eq. (11),
(42)
z \#4#6z #4z #z 8
z \#4z \#6#4z #z 8
"(z !z \) N(z ).
(43)
The 2D operator derives from these relations as it has to perform the interpolation with smoothing in one direction (relations (40) and (41)), and the interpolation with smoothing and derivation in the other direction (relations (42) and (43)). As these
A. Chehikian / Signal Processing 78 (1999) 329}347
two operations involve some common "ltering, it is possible to reduce the computation cost with the following algorithm: E Step 1. Compute the cubic spline coezcients while smoothing: y( (z , z )"s(z , z ) f (z ) f (z ). N N E Step 2. Expand the spline coezcients by inserting one column and one row of zero between each column and line: t y( (z , z ). y( (z , z ) & E Step 3. Filter the expanded coezcients in the two directions: y( (z , z )"y( (z , z )N(z )N(z ). One may notice that the sequence of steps 2 and 3, realizes the EXPAND operation as de"ned by Burt [1]. E Step 4. Compute the gradient components: s( (z , z )"K y( (z , z ) (z , z \)B(z ), N s( (z , z )"K y( (z , z ) (z , z \)B(z ). N K is a normalizing factor which ensures that the N gradient relative to a step contour of unit amplitude is equal to 1. 4(1#b) , K " (1!b)(3#b) 4(1#b) K " . (1!b)(3(1#b)#2b) For that interpolated operator, we can de"ne a detection criterion which makes more sense, in the discrete case, than the localization criterion de"ned by Canny. We consider "rst a noiseless 1D discrete step of amplitude A. The gradient estimated by our operator will exhibit a maximum of the same amplitude A located at k "k . The sharper this max imum, the less it can be disturbed by the noise.
341
A measure of the sharpness of the maximum is given by S"s( (k )!s( (k !1)"s( (k )!s( (k #1)
1!b A 3#b
if p"1,
(1!b) A 3(1#b)#2b
if p"2.
"
Sharpness depends on the smoothing parameter which also controls the noise level on the gradient. We consider now that the step is corrupted by adding a white noise of variance n. The output noise energy will be given by
>
u(x) dx, \ where u(x) is the continuous response of the proposed operator to the discrete unit impulse +i(k), as de"ned by Eq. (13). For a given value of the smoothing parameter b, we would like the sharpness measure S to be as high as possible and the output noise RMS value N as small as possible. We then consider the detection criterion D which should be as high as possible:
N"n
S S/A A A D" " "D . n N (> u(x) dx n \ This detection criterion can be split into two parts: one depending on the input: A/n , the other de pending only on the operator: D. Fig. 7 displays the variation of the intrinsic sharpness S/A and the intrinsic detection criterion D versus the signal-tonoise criterion R. It is evident from this "gure that better detection performance can be obtained with the smoothing "lter of order 1 than with the smoothing "lter of order 2. Moreover, the use of this "lter results in a lower computation cost, roughly (5 multiplies and 9 adds) times (2N;2M), where N and M are the number of rows and columns of the image. Fig. 8 shows the contours detected using the proposed "rst-order interpolated operator on an aerial image. For the sake of comparison, Fig. 9 displays the contours detected using Canny's operator
342
A. Chehikian / Signal Processing 78 (1999) 329}347
Fig. 7. Intrinsic sharpness and detection criterion versus signal-to-noise criterion. Solid line: "lter F . Dashed line: "lter F .
Fig. 8. Contours detected with the proposed operator.
(top-right), Deriche's operator (bottom-left) and Shen}Castan's operator (bottom-right). For all operators, the signal-to-noise criterion is set to the same value R"1. The amplitude and orientation of the gradient are computed as usual, and the gradient amplitude is non-maxima suppressed in the gradient direction with low and high thresholds set to 2 and 8, respectively. One can observe that even very close contours, corresponding to roads,
taxiways, gear tracks, are detected when using the operator we propose. This is mainly due to the fact that contours are detected on a double sized lattice o!ering distinct sites to contours and image pixels. Moreover, because the gradient provided by our operator exhibits a true maximum in presence of contours, the search of such maxima is easier and more e$cient than using the popular Canny's or Deriche's or Shen-Castan's operators.
A. Chehikian / Signal Processing 78 (1999) 329}347
343
Fig. 9. Contours detected with Canny's operator (top-right), Deriche's operator (bottom-left), Shen}Castan's operator (bottom-right).
5. Application to image segmentation We describe in the following, an algorithm which uses cooperation between contours and regions to get the best segmentation of an image. Because the contour detection algorithm we use gives highquality contours, image segmentation becomes trivial. The segmentation algorithm includes the following steps: E Contour detection using the previous gradient operator. This step requires selection of three control parameters. Low and high thresholds depend on the image properties and should be chosen regarding these properties. It is imperative, for the following steps of the algorithm, that the probability of failing to detect real edge points be minimized. It can be seen from Fig. 5, that increasing the signal-to-noise criterion R results in decreasing the intrinsic detection criterion D. Choosing R 1 and thus a"1.5 is a good compromise, experimentally veri"ed. With this value of the smoothing parameter, real contours are e$ciently detected but false contours due to noise are detected as well. They will be eliminated further using region properties.
E Contour closing: interrupted contours are continued along the gradient ridges until a contour is reached. Continuation is authorized inside an angle of 903 in the direction of the interrupted contour, thus avoiding the contour to turn back to the starting point. Contours may be continued in more than one direction depending on the gradient topology. Contours are closed in two steps. In the "rst step, gradient ridges are detected by looking for a local maximum of the gradient magnitude, whatever the gradient orientation. During this step, contours are prolonged as much as the gradient magnitude exhibits a ridge. This step tends to close a contour when the gradient magnitude vanishes and become lower than low threshold. In a second step, contours are prolonged for four pixels at a maximum, in the direction where the gradient magnitude exhibits the highest value. This second step tends to close a contour in the vicinity of a junction. This very simple algorithm appears to be e$cient for contour closing. E Region labeling: the interior of a closed contour is "lled up with the same label. E Contour pruning: contours that are not closed after step 2 do not de"ne a region, they can be
344
A. Chehikian / Signal Processing 78 (1999) 329}347
Fig. 10.
considered as not useful for image understanding. Such contours, having a unique region label in their 8-neighborhood, are very easily recognized. They are canceled and replaced by this region label. At the end of this step, images are usually oversegmented. This is due to the fact that image regions being not strictly uniform, the gradient magnitude does not exactly "t the contrast between the two adjacent regions. It may also be due to the contour closing algorithm. This can be corrected by a merging step: E Region merging: two adjacent regions are merged if the absolute value of their di!erence of average gray level is lower than the high threshold used for contour selection. From our experience, this algorithm gives nice results for a large variety of images, as illustrated in Figs. 10}13. In these "gures: an original image is represented on the top left; the result of segmentation is represented on the right, contours are depicted in white, while regions are depicted by their average gray level; the image reconstructed from the segmented image is represented on the bottom left, each pixel's gray level is the average gray level
of the region to which it belongs. Figs. 10 and 11 represent the segmentation results of aerial images: an airdrome (Fig. 10) and "elds (Fig. 11). Such images are usually easy to segment. Fig. 12 represents the segmentation result of the picture of a lady, while Fig. 13 represents the segmentation result of an infra-red image of buildings. These last images are more di$cult to segment: the lady image is made of non-homogenous regions and the contrast between the face and the background is very low; the infra-red image is very dark and of low contrast, it includes a lot of small regions. All images have been segmented using the same control parameters: the smoothing parameter a has been set to 1.5, low and high thresholds have been set to 2 and 8, respectively. The segmentation process takes less than 1 sec on a Sun Ultra1-166 MHz workstation, depending on the image complexity. Table 1 gives for each image, the number of regions obtained before and after the merging step along with the PPSNR of the reconstructed image. Comparing the PPSNR of the reconstructed image before and after merging, gives some cue about the choice of the low and high threshold. For `airdromea, `"eldsa and `ladya images, the merging
A. Chehikian / Signal Processing 78 (1999) 329}347
345
Fig. 11.
Fig. 12.
step reduces the number of regions by roughly 30}40%, while the PPSNR is reduced by 0.1}0.5 dB. For the `buildinga image, the number of regions is reduced by 50% and the PPSNR is
reduced by 1.7 dB, which is not negligible. As the high threshold only is used as a merging criterion, one can deduce that this control parameter is too high. Setting the high threshold to 4 results in
346
A. Chehikian / Signal Processing 78 (1999) 329}347
Fig. 13.
Fig. 14.
a number of regions before and after merging going from 3676 to 2501 (30% reduction) and a PPSNR going from 34.22 to 34.09 dB (0.13 dB reduction). Fig. 14 displays the segmentation which results.
The segmented image appears to be very complex, however it corresponds to the image content. Obviously, other a priori knowledges should be used to obtain a more signi"cant segmentation.
A. Chehikian / Signal Processing 78 (1999) 329}347
347
Table 1 Airdrome
Number of regions PPSNR
Fields
Lady
Buildings
Before merging
After merging
Before merging
After merging
Before merging
After merging
Before merging
After merging
748
524
958
579
1744
1040
2428
1245
31.7dB
31.6dB
34.0dB
6. Conclusion This paper concerns the problem of image segmentation. In the algorithm we propose, contours and regions cooperate to achieve fast and e$cient segmentation. Contours are detected "rst and regions are obtained by "lling the interior of contours by the same label. Conversely, regions are used for eliminating false contours that are generated by the noise or are insigni"cant for image understanding. From our point of view, this strategy for reducing noise e!ects is better than the usual image bandwidth reduction which always lead to information lost. Moreover, our algorithm facilitates the choice of control parameters: contrast parameters make more sense to a user than homogeneity parameters. To be e$cient, our approach requires a contour detection algorithm able to provide contour pixels located between image pixels. Therefore, we have developed a gradient operator based on B-spline interpolation, which ensures that a unique contour generate a unique response. One may also wonder if it is only a coincidence to discover with this approach, the well known Shen}Castan and Deriche operators. We also propose an interpolated version of this operator which ensures that,
33.5dB
28.1dB
27.9dB
33.4dB
31.7dB
provided that region boundaries can be modeled by step edges or have similar symmetry properties, the maximum of the gradient is exactly located midway of image pixels. Because segmentation is achieved on a double sized lattice, thus avoiding interference between image and contour pixels, contours and regions can cooperate e$ciently.
References [1] P. Burt, E.H. Adelson, The laplacian Pyramid as a compact image code, IEEE Trans. Commun COM-31 (1983) 532}540. [2] J.F. Canny, A computational approach to edge detection, IEEE Trans. Pattern Anal. Mach. Intell. PAMI-8 (1986) 679}697. [3] S. Castan, J. Zhao, J. Shen, Optimal "lter for edge detection methods and results, in: O. Faugeras (Ed.), ECCV 90, Lecture Notes in Computer Science, Springer, Berlin, 1990, pp. 13}17. [4] R. Deriche, Using Canny's criteria to derive a recursively implemented optimal edge detector, Internat. J. Comput. Vision 1 (1987) 167}187. [5] M. Unser, A. Aldroubi, M. Eden, B-Spline signal processing: Part I } Theory, IEEE Trans. Signal Process. 41 (1993) 821}833. [6] M. Unser, A. Aldroubi, M. Eden, B-Spline signal processing: Part II * E$cient design and applications, IEEE Trans. Signal Process. 41 (1993) 834}848.
Image segmentation by contours and regions cooperation A. Chehikian* Laboratoire des Images et des Signaux, Institut National Polytechnique de Grenoble and Universite& Joseph Fourier Grenoble, INPG, 46 Avenue Felix Viallet, 38031 Grenoble Cedex, France Received 13 July 1998; received in revised form 1 February 1999
Abstract In this paper we propose an algorithm for image segmentation in which contours and regions cooperate. In order for this cooperation to be e$cient, the role of contour pixels and region pixels should be considered di!erently. Region pixels carry information about the content of the image. They are de"ned on a lattice identical to the image lattice. Contour pixels carry information about regions boundary. They should be located between region pixels, therefore on a lattice distinct than the image lattice. A way to make these two statements concord is to achieve segmentation on a double-sized lattice where sites of both even coordinates are reserved for image pixels, whereas the remaining sites are reserved for contour pixels. We therefore propose a gradient operator, based on B-spline interpolation, which locates contour pixels between image pixels. Contours detected in this manner are used in an image segmentation algorithm which makes contours and region cooperate. 1999 Elsevier Science B.V. All rights reserved. Zusammenfassung In diesem Beitrag schlagen wir einen Algorithmus zur Bildsegmentierung vor, in welchem Konturen und Bereiche kooperieren. Um diese Kooperation e$zient zu gestalten, sollten die Konturpunkte und Bereichspunkte unterschiedlich behandelt werden. Bereichspunkte beinhalten Informationen uK ber den Bildinhalt. Sie sind de"niert uK ber eine Gittergleichheit zum Bildgitter. Konturpunkte beinhalten Informationen uK ber Bereichsgrenzen. Sie sollten zwischen Bereichspunkten liegen und damit auf einen vom Bildgitter unterschiedlichen Gitter. Ein Weg, diese Darstellungen in Einklang zu bringen, ist die Scha!ung der Segmentierung auf einem Gitter doppelter GroK {e, wobei PlaK tze beider geraden Koordinaten fuK r Bildpunkte und die uK brigen PlaK tze fuK r Konturpunkte reserviert sind. Wir schlagen hierzu einen Gradientenoperator basierend auf B-Spline-Interpolation vor, welcher die Konturpunkte zwischen den Bildpunkten lokalisiert. Auf diese Weise detektierte Kontouren werden in einem Bildsegmentierungsalgorithmus benutzt, welcher Konturen und Bereiche kooperieren laK sst. 1999 Elsevier Science B.V. All rights reserved. Re2 sume2 On propose dans cet article un algorithme de segmentation d'image par coopeH ration entre contours et reH gions. Pour que cette coopeH ration soit e$cace, les pixels de contour et ceux de reH gions doivent e( tre consideH reH s di!eH remment. Les pixels formant une reH gion apportent une information sur le contenu de l'image et sont donc deH "nis sur une grille d'eH chantillonnage identique a` celle de l'image. Les pixels deH "nissant les contours apportent, eux, une information sur la frontie`re des reH gions. LocaliseH s entre les reH gions, ils doivent e( tre deH "nis sur une grille d'eH chantillonnage distincte de celle
* Corresponding author. Tel.: #334-7657-4358; fax: #334-7657-4790. E-mail address: [email protected] (A. Chehikian) 0165-1684/99/$ - see front matter 1999 Elsevier Science B.V. All rights reserved. PII: S 0 1 6 5 - 1 6 8 4 ( 9 9 ) 0 0 0 7 3 - 0
330
A. Chehikian / Signal Processing 78 (1999) 329}347
de l'image. Pour assurer la coheH rence de chaque type de pixel, on propose de reH aliser la segmentation sur une grille d'eH chantillonnage de reH solution double de celle de l'image. Sur cette grille, les sites de coordonneH es toutes deux paires sont reH serveH s aux pixels d'image, les autres sont reH serveH s aux pixels de contour. Nous proposons un opeH rateur d'estimation du gradient, baseH sur l'interpolation B-spline, qui situe les contours entre les pixels d'image. Nous utilisons les contours ainsi deH tecteH s, dans un algorithme de segmentation mettant en wuvre une coopeH ration contours-reH gions. 1999 Elsevier Science B.V. All rights reserved. Keywords: Image segmentation; Contour detection; B-spline interpolation
1. Introduction Image segmentation is a necessary step in image understanding or low bit-rate coding. The goal of image segmentation is to split the image into regions. A region is a set of connected pixels having roughly the same value for some property, usually the gray level (or some other scalar property), surrounded by a closed contour representing its boundary. Two classes of image segmentation algorithms may be found in the literature. In the "rst, closed contours are detected, and regions are then deduced from their boundary by "lling the interior with the same label. Algorithms from the second class split the image into regions using some homogeneity criterion. Boundaries are then deduced if needed. Both type of algorithms have advantages and drawbacks. Contour detection is sensitive to noise and may result in both over-segmentation and under-segmentation. On the other hand, algorithms for contour detection are easy to implement because of the local nature of operations. Region based algorithms are usually more demanding in computation while being less sensitive to noise. Contour and region segmentation processes are usually considered as dual. One may take advantage of the two methods to obtain the best segmentation according to some segmentation quality criterion. However, such a cooperation requires that the role of each kind of pixel is well de"ned. Region pixels carry informations about an image property, they are therefore de"ned on the image lattice. Contour pixels carry informations about region boundaries, they should be located between regions and therefore should be de"ned on a distinct lattice.
In this paper, we propose an algorithm which makes contours and regions cooperate to obtain fast and e$cient image segmentation. In this algorithm, contours are detected "rst, regions are then obtained by "lling the interior of closed contours with the same label. The resulting regions are then used to validate signi"cant contours: close contours that separate two regions which do not satisfy some merging criterion. Final segmentation results from this cooperative process. With this approach, contour detection is a key process for which we propose an e$cient algorithm. Most contour detection algorithms are based on the evaluation of spatial gradient or Laplacian. Pixels are labeled contour pixels if an extrema of the gradient or a zero crossing of the Laplacian is observed. Such contour detection schemes raise the question of evaluating derivatives from a sampled signal, and also the question of the in#uence of noise when computing such derivatives. Canny [2], in his well-known approach, considers the case of a continuous edge modeled by a unit step. He looks for a "lter which exhibits a maximum output when such a step is present. To be optimal, such a "lter should: maximize the probability of detecting a contour, minimize the probability of detecting false contours due to noise, minimize the average displacement of the contour in the presence of noise and "nally, maximize the average distance between the detected false contours. The optimal "lter he proposes, is approximated by the "rst derivative of a Gaussian, so that contour detection reduces to searching the maximum of the "rst derivative of the output of a smoothing "lter. Deriche [4], considering an in"nite impulse response "lter, gives an exact solution for Canny's problem. Shen and Castan [3] also consider the case of a continuous step as
A. Chehikian / Signal Processing 78 (1999) 329}347
a model of an edge. They look for a "lter having the following properties: when only the step is present, the energy of the "rst derivative of the output should be maximum, while this energy should be minimum when only the noise is present. They conclude that the signal must be smoothed by a "rst-order in"nite impulse response "lter, the contour is detected when the "rst derivative of the output exhibits an extremum. It is interesting to note that all of these approaches consider the case of a continuous signal; the solutions that are proposed are extended to the discrete case by sampling the continuous impulse responses. In the contour detection algorithm we propose, we consider the more realistic case of a contour modeled by a discrete step. As we aim to estimate the "rst and/or the second derivative in the presence of such a signal, not only on sample points but between them, and because it is computationally e$cient, we have chosen to "t a continuous cubic spline to the discrete samples. However, in the presence of step contours, spline functions tend to oscillate and therefore provide multiple extrema around the true contour. It is therefore necessary to regularize the sampled signal in order to ensure the uniqueness of the gradient maximum or the Laplacian zero-crossing, while keeping a good localization of the contour. This paper is organized as follows. In Section 1, we remind some fundamental results concerning spline interpolation. These results are extracted from [5,6], readers are suggested to refer to these papers for more details. In Section 2, we consider the approximation of a discrete step by a cubic spline and deduce the values of the "rst and second derivatives at sample points. This will point out that cubic splines are not convenient for step contour detection, unless the bandwidth be reduced by an appropriate "lter. We then propose a set of regularization "lters which ensure the uniqueness of the gradient maximum and the zero-crossing of the Laplacian. We deduce the resulting 2D gradient operator, and show that a strong similarity exists between this operator and those of Shen}Castan and Deriche. However, if the continuous response underlying the discrete response of the proposed
331
operator is up-sampled on a double-sized lattice, the resulting discrete operator performs signi"cantly better than these popular edge detectors. In Section 3, we derive the interpolated operator which provides contour pixels located between the image pixels. We de"ne a detection criterion which makes more sense, in the discrete case, than the localization criterion, K, de"ned by Canny. We show that the best smoothing "lter, regarding the detection behavior, is again the one which induces the lowest computation cost. We illustrate on some examples, the bene"t we can get from our operator comparing with the well-known Deriche or Shen}Castan operators. In Section 4, "nally, we propose a very simple algorithm for image segmentation which makes contours and regions cooperate, and give some experimental results.
2. Discrete signal approximation by a continuous spline: some fundamental relations Unser et al. [5,6] consider the problem of approximating a signal discretized on a regularly spaced grid, by a spline function of degree n. Such a spline function is continuous and equal to a polynomial of degree n on each interval [k, k#1] if n is odd, on [k!, k#] if n is even (k3Z), continu ously derivable up to the order n!1, it is entirely de"ned by the set of coe$cients +y(k)3l , by the equation > gL(x)" y(k) bL(x!k), x3R, I\
(1)
where bL(x) is the continuous B-spline of order n. Such a function will exactly approximate a discrete signal +s(k), if +gL(k),"gL(x)" "+y(k),*+bL(k),"+s(k),, VI
(2)
where +bL(k), is the discrete B-spline. From Eq. (2), and considering the z-transform, we can write Eqs. (3) and (4) which de"ne a reversible transformation. Eq. (3) de"nes the direct transform giving the set of coe$cients +y(k), from the set
332
A. Chehikian / Signal Processing 78 (1999) 329}347
+s(k), whereas Eq. (4) de"nes the inverse transform giving +s(k), from the set +y(k),.
at a higher sampling rate through the use of the inverse transform:
X y(z)"s(z)BL(z)\, +y(k),&
(3)
X s(z)"y(z)BL(z). +s(k),&
(4)
+s (k),"+gL (k),"+[y(k)] ,*+bL (k),, (8) K K tK K where m is the zooming rate, +[y(k)] , is the set tK +y(k), expanded by the factor m:
In these equations, BL(z)\ is the transfer function of the direct B-spline "lter of order n, whereas BL(z) is the transfer function of the indirect B-spline "lter. A useful property of these transforms lies in the fact that the derivatives of the continuous spline being continuous splines (of lower degree) they are fully characterized by their coe$cients which can be easily deduced from the coe$cients y(k). First and second derivative values at sample points are obtained by
*gL(x) *x
and +bL (k), is the up-sampled discrete B-spline. K Eqs. (3), (5) and (6) will be used to obtain the value of the "rst and second derivative at sample points. In addition, Eq. (8) will provide the same value on a grid which is m times denser. To use these relations we need to know the coe$cients of the direct and indirect "lters. They are obtained from the recursive equations:
(10)
((n#2)/2#k/m)bL\(k#m)#(n/2!k/m)bL\(k) K K , cL (k)" K n
(11)
(5)
"+y(k#1)!2y(k) VI #y(k!1),*+bL\(k),,
with the starting conditions: m m 1 for ! )k) , 2 2 b (k)" K 0 else,
(6)
1 for 1!m)k)0, c (k)" K 0 else.
(7)
3. Edge detection using the cubic spline transform
where +cL\(k), is the discrete shifted B-spline: cL\(k)"bL\(k#).
(9)
((n#1)/2#k/m)cL\(k)#((n#1)/2!k/m)cL\(k!m) K K bL (k)" , K n
"+y(k)!y(k!1),*+cL\(k),, VI
*gL(x) *x
y(k) if k "mk, +[y(k)] ,"y(k )" tK 0 elsewhere,
These relations should be interpreted as follows. The coe$cients of the "rst derivative are obtained by the di!erence [y(k)!y(k!1)], values of the "rst derivative at sample points are further obtained from the convolution with the shifted Bspline of degree n!1. In the same way, the coe$cients of the second derivative are obtained by [y(k#1)!2y(k)#y(k!1)], samples of the second derivative are then obtained by convolving with the symmetric B-spline of degree n!2. Another useful property of the spline transform lies in the fact that the signal may be reconstructed
3.1. One-dimensional formulation In the following, we consider a discrete 1D step contour modeled by a discrete unit step +u(k),:
u(k)"
1 for k*0,
(12)
0 else.
We shall also consider the unit impulse +i(k),:
i(k)"u(k)!u(k!1)"
1 if k"0, 0 else.
(13)
A. Chehikian / Signal Processing 78 (1999) 329}347
As we are only interested in the "rst and second derivative at sample points, a cubic spline transform is quite appropriate and results in lower computation cost than if a higher degree spline transforms were used. In order to compute the values of the "rst and second derivatives, we need to know the transfer functions of the direct spline "lter of order 3 to obtain the spline coe$cients, the indirect shifted "lter of order 2 to obtain the sampled "rst derivative, and the indirect "lter of order 1 to obtain the sampled second derivative. Parameters of these transfer function are obtained from Eqs. (10) and (11): B(z)"1,
z#1 , C(z)" 2
One can notice that if indirect "lters are of "nite response, the direct spline "lter is of in"nite impulse response. Writing 6 B(z)\" z#4#z\
1 1!a 1 " # !1 , 1#a 1!az 1!az\ where a"(3!2(0, it is easy to show that the discrete impulse response +r(k), of the direct spline "lter of order 3 is given by 1!a r(k)" a I . 1#a
(14)
Eq. (3) yields +y(k),"+u(k),*+r(k),.
(15)
It results from Eqs. (5) and (6) that the spline coe$cients of the "rst and second derivatives are given by Eqs. (16) and (17), respectively: +y(k)!y(k!1),"+u(k)!u(k!1),*+r(k), "+i(k),*+r(k),"+r(k),,
+y(k#1)!2y(k)#y(k!1),"+r(k#1)!r(k),. (17) Therefore, the values of the "rst and second derivatives at sample points are given by
+u (k),"
r(k)#r(k#1) 2
(18)
(19)
1!a aI for k*0, 2 Nu (k)" 1!a a\I> for k(0, 2 +u (k),"+r(k#1)!r(k),
6 B(z)\" . z#4#z\
333
(16)
(1!a) ! aI for k*0, (1#a) Nu (k)" (1!a) a\I> for k(0. # (1#a)
Fig. 1(a) displays the continuous cubic spline approximating the discrete 1D step, as well as its "rst derivative, whereas Fig. 1(b) displays the second derivative. Sample points are marked with `a. We can see on these "gures that, in the vicinity of the step, the maximum of the "rst derivative and the zero crossing of the second derivative occur exactly at the midpoint of the gray-level transition, i.e. between the two image pixels. However, these derivatives exhibit multiple extrema and zero-crossings, respectively. Cubic B-spline interpolation cannot be used in this way for contour detection unless the signal bandwidth be reduced by an appropriate smoothing "lter.
3.2. Image bandwith reduction for contour detection by cubic spline transform We need to de"ne a smoothing "lter of impulse response +F (k),, which ensures that the second N derivative of the cubic spline which interpolates the discrete step smoothed by this "lter, will exhibit a unique zero crossing. Looking for a "lter of
334
A. Chehikian / Signal Processing 78 (1999) 329}347
Fig. 1. (a) First and (b) second derivative of the cubic spline approximating a step edge.
E The zero crossing of the Laplacian should be, for symmetry reasons, located at x"!0.5:
transfer function of the form X F (z)"B(z) f (z) +F (k),& N N N
u( (!1)"!u( (0)Nf (1)"f (!1). N N
will simplify the problem. Indeed, the coe$cients of the cubic spline interpolating the smoothed signal, will be given by X y( (z)"s(z)F (z)B(z)\"s(z) f (z) +y( (k),& N N
E The zero crossing of the Laplacian should be unique: u( (k))0, ∀k*0Nf (k#1)!f (k))0, ∀k*0, N N u( (k)*0, ∀k(0Nf (k#1)!f (k)*0, ∀k(0. N N (22)
X +s(k), +f (k), & * N and the values of the "rst and second derivative of the smoothed step +u( (k), at sample points will be, from Eqs. (18) and (19): f (k#1)#f (k) N , u( (k)" N 2
These three conditions are satis"ed by any "lter of positive, symmetric, unimodal impulse response: a Gaussian "lter for example, would be a good candidate. However, because of their interesting properties, we focus on the family of "lters of transfer function:
N 1!b 1!b f (z)" , N 1!bz 1!bz\
u( (k)"f (k#1)!f (k). N N
b"e\?, a'0Nb3]0, 1[ ,
We can now express some conditions that are to be satis"ed.
whose impulse response are of the form
E The gradient estimated on a positive step contour, should be positive:
+ f (k)," N
f (k)*0, N
∀k.
(21)
(20)
1!b N "k"N\ N\ b I # c "k"G G 1#b (p!1)! G
c '0 ∀i. G
(23)
,
(24)
A. Chehikian / Signal Processing 78 (1999) 329}347
Such impulse responses are positive, symmetric, and unimodal and therefore satisfy conditions (20)}(22). The "rst two members of this family of "lters, as it will be seen in the following, lead to a very popular contour detection operators. Their transfer function and impulse response are given by 1!b 1!b f (z)" , 1!bz 1!bz\
1!b b I , + f (k)," 1#b
(25)
(26)
335
1!2b#b 1!2b#b , f (z)" 1!2bz\#bz\ 1!2bz#bz
+ f (k),"
1#b 1!b b I "k"# 1#b 1!b
.
(27)
(28)
Fig. 2 displays the continuous cubic spline which approximates the discrete step smoothed by the "lter F , along with its "rst derivative (Fig. 2(a)) and its second derivative (Fig. 2(b)). Fig. 3 displays the same results when using the "lter F instead of F . In both "gures, a is set to 2. One can observe from these "gures that the two "lters, even of di!erent order, provide continuous
Fig. 2. (a) First and (b) second derivative of the cubic spline approximating a step edge smoothed by F .
Fig. 3. (a) First and (b) second derivative of the cubic spline approximating a step edge smoothed by F .
336
A. Chehikian / Signal Processing 78 (1999) 329}347
splines and spline derivatives of the same shape: piecewise cubic for the spline, piecewise quadratic for its "rst derivative and piecewise linear for its second derivative. One may wonder if the use of a high-order smoothing "lter, which of course induces a higher computation cost, is more bene"cial than the use of a lower-order smoothing "lter. This point will be discussed later. We would like to "rst de"ne the 2D operators resulting from the use of these smoothing "lters. 3.3. Deriving a 2D gradient operator When using the smoothing "lters de"ned in the previous section, it results from relations (3), (5) and (4) that the gradient values at sample points are obtained by
whereas the transfer functions of the Shen}Castan "lters may be expressed by
1!b 1!b , F (z)" 1 1!bz 1!bz\
One can see that the derivation "lters di!er only by a normalizing factor, whereas the smoothing "lters di!er by the indirect spline "lter which is not present in the Shen}Castan's operator. Using the smoothing "lter of order 2: z#4#z\ F (z)" 6 ;
(36)
;
where (z!z\) H (z)"f (z) N N 2
s( (z , z )"s(z , z )H (z )F (z ), (30) N N s( (z , z )"s(z , z )H (z )F (z ). (31) N N These relations are very similar to those proposed by Shen}Castan and Deriche for the estimation of the gradient vector components. It is worthwhile to compare the proposed operator to theirs. Using the smoothing "lter of order 1:
z#4#z\ 1!b 1!b F (z)" , 6 1!bz 1!bz\
1!2b#b 1!2b#b , 1!2bz#bz 1!2bz\#bz\ (37)
(29)
and p is the order of the smoothing "lter. In the 2D case, smoothing is performed in the two orthogonal directions, whereas derivation is performed only in the direction of the gradient's component. It results that the two gradient components will be obtained by
1!2b#b 1!2b#b , 1!2bz#bz 1!2bz\#bz\
(z!z\) H (z)" 2
(z!z\) s( (z)"s(z) f (z) "s(z) H (z), N N 2
(z!z\) 1!b 1!b H (z)" , 2 1!bz 1!bz\
b(1#b) 1!b 1!b H (z)" (z!z\) . (35) 1 (1!b) 1!bz 1!bz\
X s( (z)"[[s(z) F (z)]B(z)\] (1!z\) (1#z), s( (k)& N 2
(34)
whereas the transfer functions of the Deriche "lters (case u equal to 0) may be expressed by F (z)"(c z#c #c z\) " 1!2b#b 1!2b#b ; , 1!2bz#bz 1!2bz\#bz\
(38) 1#b H (z)" (z!z\) " 1!b ;
1!2b#b 1!2b#b , 1!2bz#bz 1!2bz\#bz\ (39)
where
(32)
b (b\!4a!b) c " , (1!b) (b\#2a!b)
(33)
b ((a!1)b\#(a#1)b) c " . (1!b) (b\#2a!b)
A. Chehikian / Signal Processing 78 (1999) 329}347
Once more, the derivation "lters di!er only by a normalizing factor, whereas the smoothing "lters are very similar. Moreover, it can be shown that when a goes to zero, the coe$cients c and c in the Deriche smoothing "lter go to 4/6 and 1/6, respectively (a equal to 1 yields c and c equal to 0.6886 and 0.1557, respectively). One may conclude that our proposed operator and Deriche's operator, although developed from two very di!erent approaches, are identical when smoothing becomes e!ective. However, we must repeat that these comparisons make sense only when considering the discrete "lters, as the continuous "lters are quite di!erent. We compare in the next section, the response of Deriche and Shen}Castan "lters to those of the operators we propose.
3.4. Comparing with Deriche and Shen}Castan xlters A common characteristic of these "lters is that they are de"ned to optimize some criterions when considering a contour modeled by a continuous step. The results which are obtained are then extended to the discrete case by sampling the continuous impulse response. Figs. 4 and 5 show the response of the continuous and discrete "lters when applied to a step contour without and with an additive noise. The smoothing parameter of each "lter has been chosen to ensure the same signal-tonoise criterion R"1, as de"ned in Canny's approach. In both "gures, the ratio A/n of the step amplitude to the noise standard deviation is 8. One may observe in Fig. 4 illustrating the continuous case that, for both "lters, even in the presence of noise, the maximum of the response is quite well localized, due to the moderate noise environment. Considering the discrete case, illustrated in Fig. 5, results in two consequences. First, due to the violation of the sampling theorem, the exact position of the contour is lost. Indeed, any continuous step localized at x 3]!1, 0] will be described by the same discrete step. Second, the response of the discrete "lters exhibits a maximum in form of a plateau, when applied to an un-noisy step, whereas it exhibits a true maximum, but randomly
337
located at x"0 or !1, when applied to a moderately noisy step. As a matter of fact, extending Canny's localization criterion, initially de"ned from a continuous model of contour, to the discrete case is at least improper. Fig. 6 shows the normalized response of our proposed "rst-order "lter. Because in our approach the contour we consider is modeled by a discrete step, we only consider the response of the "lter to such a discrete step. The sets of points marked with `a constitute the response of the "lter discretized on the input lattice. When considering these responses, one may conclude, in accordance with Eqs. (33) and (35), that the proposed "lter exhibits the same behavior as the one of Shen}Castan (the same conclusion would be derived from a comparison between our second-order "lter and the one of Deriche). However, because we interpolate the discrete step by a continuous spline, we are able to provide the underlying continuous response and the up-sampled discrete response as well. The sets of points marked with ` 6 *a constitute the discrete response up-sampled on a double sized lattice. Such response exhibit a true maximum at x"!0.5, the position of which is very less a!ected by the noise than the response of the "lter discretized on the input lattice. The following section de"nes an interpolated operator for contour detection which reduces the uncertainty on contour position and provides contour sites located between the image pixels, thus facilitating image segmentation.
4. An interpolated operator for contour detection The structure of such an operator is very similar to the one of the previous 2D operator except that it includes an interpolating step according to Eqs. (8) and (9). For the sake of simplicity, we de"ne "rst the operations to perform to interpolate without derivation and next when interpolating with the "rst-order derivation. From Eq. (8), we can write +s( (k),"+[y( (k)] ,*+b(k), t t X ' & s( (z )"y( (z ) B(z ), t
(40)
338
A. Chehikian / Signal Processing 78 (1999) 329}347
Fig. 4. Shen}Castan and Deriche "lters response to a continuous step.
A. Chehikian / Signal Processing 78 (1999) 329}347
Fig. 5. Shen}Castan and Deriche "lters response to a discrete step.
339
340
A. Chehikian / Signal Processing 78 (1999) 329}347
Fig. 6. The proposed "rst order "lter's response to a discrete step.
with from Eq. (10), 0.125z \#z\#2.875z \#4#2.875z #z #0.125z B(z )" 6 which can be factorized into z \#4#z z\#4z \#6#4z #z B(z )" 6 8
(1!z \) C(z )"(1!z \) can be factorized into
(41)
(1!z \)C(z )
where one may recognize (up to a multiplicative factor 2) in N(z ) the binomial "lter of order 4. In the same manner, from Eqs. (5) and (8):
"(z !z \)
"B(z ) N(z ),
+s( (k),"+[y( (k)!y( (k!1)] ,*+c(k),, t t s( (z )"y( (z ) (1!z \) c(z ), t where, from Eq. (11),
(42)
z \#4#6z #4z #z 8
z \#4z \#6#4z #z 8
"(z !z \) N(z ).
(43)
The 2D operator derives from these relations as it has to perform the interpolation with smoothing in one direction (relations (40) and (41)), and the interpolation with smoothing and derivation in the other direction (relations (42) and (43)). As these
A. Chehikian / Signal Processing 78 (1999) 329}347
two operations involve some common "ltering, it is possible to reduce the computation cost with the following algorithm: E Step 1. Compute the cubic spline coezcients while smoothing: y( (z , z )"s(z , z ) f (z ) f (z ). N N E Step 2. Expand the spline coezcients by inserting one column and one row of zero between each column and line: t y( (z , z ). y( (z , z ) & E Step 3. Filter the expanded coezcients in the two directions: y( (z , z )"y( (z , z )N(z )N(z ). One may notice that the sequence of steps 2 and 3, realizes the EXPAND operation as de"ned by Burt [1]. E Step 4. Compute the gradient components: s( (z , z )"K y( (z , z ) (z , z \)B(z ), N s( (z , z )"K y( (z , z ) (z , z \)B(z ). N K is a normalizing factor which ensures that the N gradient relative to a step contour of unit amplitude is equal to 1. 4(1#b) , K " (1!b)(3#b) 4(1#b) K " . (1!b)(3(1#b)#2b) For that interpolated operator, we can de"ne a detection criterion which makes more sense, in the discrete case, than the localization criterion de"ned by Canny. We consider "rst a noiseless 1D discrete step of amplitude A. The gradient estimated by our operator will exhibit a maximum of the same amplitude A located at k "k . The sharper this max imum, the less it can be disturbed by the noise.
341
A measure of the sharpness of the maximum is given by S"s( (k )!s( (k !1)"s( (k )!s( (k #1)
1!b A 3#b
if p"1,
(1!b) A 3(1#b)#2b
if p"2.
"
Sharpness depends on the smoothing parameter which also controls the noise level on the gradient. We consider now that the step is corrupted by adding a white noise of variance n. The output noise energy will be given by
>
u(x) dx, \ where u(x) is the continuous response of the proposed operator to the discrete unit impulse +i(k), as de"ned by Eq. (13). For a given value of the smoothing parameter b, we would like the sharpness measure S to be as high as possible and the output noise RMS value N as small as possible. We then consider the detection criterion D which should be as high as possible:
N"n
S S/A A A D" " "D . n N (> u(x) dx n \ This detection criterion can be split into two parts: one depending on the input: A/n , the other de pending only on the operator: D. Fig. 7 displays the variation of the intrinsic sharpness S/A and the intrinsic detection criterion D versus the signal-tonoise criterion R. It is evident from this "gure that better detection performance can be obtained with the smoothing "lter of order 1 than with the smoothing "lter of order 2. Moreover, the use of this "lter results in a lower computation cost, roughly (5 multiplies and 9 adds) times (2N;2M), where N and M are the number of rows and columns of the image. Fig. 8 shows the contours detected using the proposed "rst-order interpolated operator on an aerial image. For the sake of comparison, Fig. 9 displays the contours detected using Canny's operator
342
A. Chehikian / Signal Processing 78 (1999) 329}347
Fig. 7. Intrinsic sharpness and detection criterion versus signal-to-noise criterion. Solid line: "lter F . Dashed line: "lter F .
Fig. 8. Contours detected with the proposed operator.
(top-right), Deriche's operator (bottom-left) and Shen}Castan's operator (bottom-right). For all operators, the signal-to-noise criterion is set to the same value R"1. The amplitude and orientation of the gradient are computed as usual, and the gradient amplitude is non-maxima suppressed in the gradient direction with low and high thresholds set to 2 and 8, respectively. One can observe that even very close contours, corresponding to roads,
taxiways, gear tracks, are detected when using the operator we propose. This is mainly due to the fact that contours are detected on a double sized lattice o!ering distinct sites to contours and image pixels. Moreover, because the gradient provided by our operator exhibits a true maximum in presence of contours, the search of such maxima is easier and more e$cient than using the popular Canny's or Deriche's or Shen-Castan's operators.
A. Chehikian / Signal Processing 78 (1999) 329}347
343
Fig. 9. Contours detected with Canny's operator (top-right), Deriche's operator (bottom-left), Shen}Castan's operator (bottom-right).
5. Application to image segmentation We describe in the following, an algorithm which uses cooperation between contours and regions to get the best segmentation of an image. Because the contour detection algorithm we use gives highquality contours, image segmentation becomes trivial. The segmentation algorithm includes the following steps: E Contour detection using the previous gradient operator. This step requires selection of three control parameters. Low and high thresholds depend on the image properties and should be chosen regarding these properties. It is imperative, for the following steps of the algorithm, that the probability of failing to detect real edge points be minimized. It can be seen from Fig. 5, that increasing the signal-to-noise criterion R results in decreasing the intrinsic detection criterion D. Choosing R 1 and thus a"1.5 is a good compromise, experimentally veri"ed. With this value of the smoothing parameter, real contours are e$ciently detected but false contours due to noise are detected as well. They will be eliminated further using region properties.
E Contour closing: interrupted contours are continued along the gradient ridges until a contour is reached. Continuation is authorized inside an angle of 903 in the direction of the interrupted contour, thus avoiding the contour to turn back to the starting point. Contours may be continued in more than one direction depending on the gradient topology. Contours are closed in two steps. In the "rst step, gradient ridges are detected by looking for a local maximum of the gradient magnitude, whatever the gradient orientation. During this step, contours are prolonged as much as the gradient magnitude exhibits a ridge. This step tends to close a contour when the gradient magnitude vanishes and become lower than low threshold. In a second step, contours are prolonged for four pixels at a maximum, in the direction where the gradient magnitude exhibits the highest value. This second step tends to close a contour in the vicinity of a junction. This very simple algorithm appears to be e$cient for contour closing. E Region labeling: the interior of a closed contour is "lled up with the same label. E Contour pruning: contours that are not closed after step 2 do not de"ne a region, they can be
344
A. Chehikian / Signal Processing 78 (1999) 329}347
Fig. 10.
considered as not useful for image understanding. Such contours, having a unique region label in their 8-neighborhood, are very easily recognized. They are canceled and replaced by this region label. At the end of this step, images are usually oversegmented. This is due to the fact that image regions being not strictly uniform, the gradient magnitude does not exactly "t the contrast between the two adjacent regions. It may also be due to the contour closing algorithm. This can be corrected by a merging step: E Region merging: two adjacent regions are merged if the absolute value of their di!erence of average gray level is lower than the high threshold used for contour selection. From our experience, this algorithm gives nice results for a large variety of images, as illustrated in Figs. 10}13. In these "gures: an original image is represented on the top left; the result of segmentation is represented on the right, contours are depicted in white, while regions are depicted by their average gray level; the image reconstructed from the segmented image is represented on the bottom left, each pixel's gray level is the average gray level
of the region to which it belongs. Figs. 10 and 11 represent the segmentation results of aerial images: an airdrome (Fig. 10) and "elds (Fig. 11). Such images are usually easy to segment. Fig. 12 represents the segmentation result of the picture of a lady, while Fig. 13 represents the segmentation result of an infra-red image of buildings. These last images are more di$cult to segment: the lady image is made of non-homogenous regions and the contrast between the face and the background is very low; the infra-red image is very dark and of low contrast, it includes a lot of small regions. All images have been segmented using the same control parameters: the smoothing parameter a has been set to 1.5, low and high thresholds have been set to 2 and 8, respectively. The segmentation process takes less than 1 sec on a Sun Ultra1-166 MHz workstation, depending on the image complexity. Table 1 gives for each image, the number of regions obtained before and after the merging step along with the PPSNR of the reconstructed image. Comparing the PPSNR of the reconstructed image before and after merging, gives some cue about the choice of the low and high threshold. For `airdromea, `"eldsa and `ladya images, the merging
A. Chehikian / Signal Processing 78 (1999) 329}347
345
Fig. 11.
Fig. 12.
step reduces the number of regions by roughly 30}40%, while the PPSNR is reduced by 0.1}0.5 dB. For the `buildinga image, the number of regions is reduced by 50% and the PPSNR is
reduced by 1.7 dB, which is not negligible. As the high threshold only is used as a merging criterion, one can deduce that this control parameter is too high. Setting the high threshold to 4 results in
346
A. Chehikian / Signal Processing 78 (1999) 329}347
Fig. 13.
Fig. 14.
a number of regions before and after merging going from 3676 to 2501 (30% reduction) and a PPSNR going from 34.22 to 34.09 dB (0.13 dB reduction). Fig. 14 displays the segmentation which results.
The segmented image appears to be very complex, however it corresponds to the image content. Obviously, other a priori knowledges should be used to obtain a more signi"cant segmentation.
A. Chehikian / Signal Processing 78 (1999) 329}347
347
Table 1 Airdrome
Number of regions PPSNR
Fields
Lady
Buildings
Before merging
After merging
Before merging
After merging
Before merging
After merging
Before merging
After merging
748
524
958
579
1744
1040
2428
1245
31.7dB
31.6dB
34.0dB
6. Conclusion This paper concerns the problem of image segmentation. In the algorithm we propose, contours and regions cooperate to achieve fast and e$cient segmentation. Contours are detected "rst and regions are obtained by "lling the interior of contours by the same label. Conversely, regions are used for eliminating false contours that are generated by the noise or are insigni"cant for image understanding. From our point of view, this strategy for reducing noise e!ects is better than the usual image bandwidth reduction which always lead to information lost. Moreover, our algorithm facilitates the choice of control parameters: contrast parameters make more sense to a user than homogeneity parameters. To be e$cient, our approach requires a contour detection algorithm able to provide contour pixels located between image pixels. Therefore, we have developed a gradient operator based on B-spline interpolation, which ensures that a unique contour generate a unique response. One may also wonder if it is only a coincidence to discover with this approach, the well known Shen}Castan and Deriche operators. We also propose an interpolated version of this operator which ensures that,
33.5dB
28.1dB
27.9dB
33.4dB
31.7dB
provided that region boundaries can be modeled by step edges or have similar symmetry properties, the maximum of the gradient is exactly located midway of image pixels. Because segmentation is achieved on a double sized lattice, thus avoiding interference between image and contour pixels, contours and regions can cooperate e$ciently.
References [1] P. Burt, E.H. Adelson, The laplacian Pyramid as a compact image code, IEEE Trans. Commun COM-31 (1983) 532}540. [2] J.F. Canny, A computational approach to edge detection, IEEE Trans. Pattern Anal. Mach. Intell. PAMI-8 (1986) 679}697. [3] S. Castan, J. Zhao, J. Shen, Optimal "lter for edge detection methods and results, in: O. Faugeras (Ed.), ECCV 90, Lecture Notes in Computer Science, Springer, Berlin, 1990, pp. 13}17. [4] R. Deriche, Using Canny's criteria to derive a recursively implemented optimal edge detector, Internat. J. Comput. Vision 1 (1987) 167}187. [5] M. Unser, A. Aldroubi, M. Eden, B-Spline signal processing: Part I } Theory, IEEE Trans. Signal Process. 41 (1993) 821}833. [6] M. Unser, A. Aldroubi, M. Eden, B-Spline signal processing: Part II * E$cient design and applications, IEEE Trans. Signal Process. 41 (1993) 834}848.