Unsupervised region-based image segmentation using texture statistics and level-set methods
I. Karoui , R. Fablet , J.M. Boucher , J.M. Augustin
GET, ENST Bretagne, CNRS TAMCIC, CS 83818 - 29238 Brest Cedex, France Ifremer, TSI-STH, Technopˆole Brest Iroise, 29280 PLOUZANE, France
Abstract – We propose a novel unsupervised region based criterion for multi-class texture segmentation. The proposed criterion relies on the maximization of a weighted sum of Kullback-Leibler measure between distributions of local texture features associated to the different image regions. Hence, the segmentation issue is stated as the maximization of the proposed criterion and a regularization term that imposes smoothness and regularity of region boundaries. The proposed approach is based on curve evolution techniques and is implemented using level-set methods. Curve evolution equations are expressed using shape derivative tools. As an application, we have tested the method using cooccurrence distributions, distributions of Gabor filter responses and wavelet packet to segment synthetic mosaics of textures from the Brodatz album, as well as real textured sonar images. Keywords – Active regions, level set, texture, unsupervised segmentation.
I. I NTRODUCTION Texture segmentation has long been an important topic in image processing. It aims at segmenting a textured image into several regions with the same texture features. An effective and efficient texture segmentation method will be very useful in applications like the analysis of aerial, biomedical and seismic images as well as the automation of industrial applications. Recently, features computed as statistics of local filter responses have been largely used in supervised texture analysis and segmentation and several studies have shown the relevance of marginals of a large set of filters to characterize textures [1], [2], [3], [4], [5]. Zhu et al. [6] proposed a maximum entropy theory for learning probabilistic texture models from a set of empirical distribution of filter responses. Gimel’farb used the difference co-occurrence statistics to model texture [7] and later, Xiuwen Liu et al. [1] proposed a local spectral histogram, defined as the marginal distributions of feature statistics for texture classification. In a supervised context, these features extracted from known textures are compared to features extracted from an image sample to be classified by a weighted Kullback-Leibler
measure [8]. Unsupervised segmentation does not require training and prior knowledge. It is generally based on the optimization of a Bayesian criterion [9], [10] or an information-theoretic based criterion within variational inference [11], [12], [13], [14]. Variational framework, especially active contour based approaches, are more appropriate for texture segmentation because they offer an efficient manner to cope with texture and geometrical features at the region level. Existing curve evolution based unsupervised segmentation methods are mainly based on the maximization of the mutual information between the region labels and the image pixel intensities [11], the maximization of the entropy associated to grey value statistics [12] or Minimum Description Length (MDL) criterion [13], [14]. In this work, we treat unsupervised multi-class texture segmentation and we aim to combine the use of non parametric marginal distributions of filter responses for texture characterization and a novel region-based criterion for image segmentation within active contour approach associated with a level-set setting. The criterion consists in maximizing the separation between different region textural statistics according to the Kullback-Leibler divergence [15]. Our method is different from the existing segmentation methods based on active contours in two major ways. First, unlike existing methods based on mutual information maximization, entropy maximization or MDL criterion, our approach is region based and in opposite to those methods, does not assume independent and identically distributed (i.i.d.) variables. Second our method is not restricted to first order grey value features: we take into account different texture features computed as marginals of texture responses to a wide set of filter and we include an automatic feature selection step. The presented paper is organized as follows. In Section 2, texture characterization and modeling are described. Section 3 briefly reviews active contour methods and the level set approach. In Section 4, the energy functional and its derivation are presented and the curve evolution equations are computed.
Section 5 shows some results on synthetic images containing textures from the Brodatz album and on natural sonar images.
that drives the active contour to a minimum of a functional. The PDE is generally derived from an energy criterion as follows: W
II. T EXTURE
Features computed as statistics of local filter responses are widely used for supervised texture analysis and segmentation [1], [2], [3], [4], [5]. Many statistical and filtering approaches have been compared [16]. Among the most effective features are co-occurrence matrices, wavelet frames, quadrature mirror filterbanks (QMF) and Gabor filters. It should be noted that none of these feature classes outperforms the others for all textures. Here, we choose a set of different filters computed for different parameters. Let , be the image of filter response to the filter taking its values in . is a set of , is the filter output value dimension and let , be their respective distributions. Using Parzen "!#%$ for the distribution window estimation [17], is computed on a domain & as: &
!
)
&
- /21
*,+.-0/21 )
T9'\[9!
a
(1)
=@>
' ! ?
A3 !CBD ?
-FE
A3d'VUV!fe][ 3g$ & if g A3d'VUV! [ 3 $ ` if _ A3d'VUV!fh][ `
bc `
(5)
otherwise e6[
34!#56#!8793;:
is a Gaussian kernel of mean 0 and variance < . where To compare texture features, we use the Kullback-Leibler divergence [15] which is a relevant similarity measure between distributions. The Kullback-Leibler distance (=%> ) between two distributions ? and , is defined as: *
(4)
]S(^ an initial curve defined by the user. is the with S Z velocity vector and X is the inward normal of S . T9'0UV! The parametric representation of the curve S is unsuitable for many applications since it does not allow for automatic change of topology, such as merging and breaking. Level set methods, introduced in Osher and Sethian [19] to track moving interfaces in the community of fluid dynamics, circumvent these topological problems. The key idea of level W set methods is to represent the evolving curve _@ & with an implicit Lipschitz function ` which is defined by:
('
T9'VUV! Z W U YX S
DESCRIPTION
A34! 793 ? A3 !HG
The region & is entirely described by the level sets ` and its geometrical quantities can be expressed by ` . Evolving the curve S in its normal direction with speed
amounts to solving the differential equation [19]: W
?
I ,
) `
'
`
[8'3;'j!
k`"^
3;'j!
(6)
where ` ^ is the initial contour.
(2)
I and ?J For two distribution sets we define their dissimilarity measure as follows:
)i
` W % U
IV. S EGMENTATION METHOD DESCRIPTION The segmentation issue is stated as the minimization of an energy criterion lml Nn l Nn lpo , where l is a texture-based data-driven term, l a regularization term and l o a term needed to cope with multi-class segmentation.
=@>(K
?
' !
L 8' 0! ,M =@> ?
(3) A. Functional terms
L
l is the data-term, it is expressed as the distance between the statistics estimated on the different image regions. The distance is evaluated according to the proposed texture similarity measure =@>qK (equation (3))
are weights such that O . The M NM weights are computed according to the discrimination powers of the associated features. In previous work [18] we proposed a Kullback-Leibler margin maximization criterion for weight estimation in supervised case. Here, we generalize the criterion to unsupervised case by the alternation between weight estimation and image segmentation. Initially the weights are set to QP9R and are readjusted during the M segmentation process. where
l
! & r r s
L
5
=@>(K t rQu rwvyxzu rwv\ { r
s
& r
!|'
& rwv
!!
(7)
!
& r is a set of filter response distributions estimated where globally inside the region & r . Because of the convexity of the Kullback-Leibler divergence, the maximization of l leads to the partition that gives homogeneous regions [20]. l penalizes the region contour lengths and is expressed by:
III. ACTIVE
CONTOURS AND
L EVEL SET IMPLEMENTATION
The idea behind active contour segmentation methods is toT T9'UV! evolve a parametric curveU S in the image domain & . may be its arc-length and is an evolution parameter. The curve evolution is described by a partial differential equation (PDE)
l
L
s
r r ,}
)
_ r
) ' }
r
$
(8)
where _ r denotes the contour of the region & r . Using level set functions and regularized Heaviside and delta functions, l can be written as follows (see [21] for details):
l
L
s
} r
r
7
* +
! )i
`#r
) 703
` r
a
b
c
E
3
n
) 3 )
if
(10)
l o
` r
!#5
793
#"
(12)
'$
The evolution equation of ` r global functional l is given by: W
` r W U
l
W n
l
W n
l o
': : :y'
'&%' $
=
': : :y'
related to the =
(13)
W
- /Q1 G98
W
& r v
!V!
TQ!V!/:;
7
)i
6 L 5 ! ' V! ! ) = >(K @ & r & r v & r r v{ r ! ! E - - /Q1 5 TQ!V! :; BD & r v n & r f p ! ! G 8 & r v & r M )
L ` r
)
(17) C. Initialization We consider a fully unsupervised case; neither the class number texture feature prototypes are known. Let ,=?>! nor Q the N,=?>|! >+@0A u be the set of features estimated T locally within a window centered at each image pixel . To initialize the segmentation map and to determine the class number, we use a k-means algorithm with an entropy based prior proposed in [9] and generalized in [10] for fuzzy k-means. The difference here with those algorithms consists in the use of the similarity measure =%> K (see equation (3)) instead of the Euclidean distance used in [9], [10]. The objective function associated to the initialization algorithm is the following: B
' '+*
` r
! '
L
s L
>+@0A =@>(K r C
& r
! ' -=D>|!!#5?EF r
BD
-
GF r
!-H
(18)
is', the evolution equation term related to the ' '-* functional term l ( . W ' '/* l.( are directly estimated from level set functions [25], [21].
& r
B. Computation of the evolution equation
l
! & r v n ! & r
E
` r W U n
r
(11)
if
* +! L s
=%> K
where X r is the unit inward normal to p & r , < X its area element W and 8 is the convolution operator. l is then given by:
segmentation in order to fulfill the image partition condition. In the literature, there are several techniques for dealing with the representation of the different classes and their boundaries by level-sets [22], [23], [24]. Our multi class model is inspired by the work of C. Zhao et 'wal. [22]. The idea is to associate a level' : :y: ' = to each region & r and the image set function ` r partition condition is expressed by the following term:
W
& r
)
W
3 ) 3 ) n DHT
if )3 ) h [
Y
*35476
Z
l o is an additional term, required to cope with multi-class
'
)
(16)
3 n T G
3 e if [ 3 h 5p
if
l)(
L
X
A3 !
W
' 2
& r
in the direction X is given by the following
r v{ r ! 5 5 L B D & r ! M & r v 2 Z 7 < TQ! X Xr X
(9)
and are Heaviside and delta so functions respectively, [ v that: when , and , and (in the sense of distributions).
W
l
A34!
2
derivative of l equation:
!
! l o ` r
r }
5
!870 1
` r
` r
` r
)i L
!
i E
r
s
` r
)G
(14)
F
is the class number, r
2P
) )L
tGIKJ$L M IKN u rwx and E a >+@0A
positive constant. For this initialization algorithm, the weights are set to 2P9R . The clustering algorithm is given below: 1. Initialize = and each class feature centroids ! & r r I s ; 2. Update the segmentation map according to the following classification rule:
,
` r
!#5
! "
(15)
The evolution equation related to l is more complex, since it involves computations over the spatial support of each region. To differentiate l , we use shape derivative tools, especially the Gˆateaux derivative theorem given in [26]. The Gˆateaux
where =
B S0TVUXWZY\[ =@>(K r
3. Update texture centroids current segmentation map;
&pr
& r
! ' ,=?>|!!
! r I s
(19)
according to the
50
50
100
100
150
150
200
200
250
250 50
100
150
200
250
Fig. 1. 3 B RODATZ TEXTURE MOSAIC UNSUPERVISED SEGMENTATION .
50
100
150
200
250
Fig. 2. 5 B RODATZ TEXTURE MOSAIC UNSUPERVISED SEGMENTATION .
F
4. Update the prior probability r of each class and discard F h the class if r ; n and return to step until convergence. 5. = m= [ The weights are readjusted every iterations of the application of the PDEs given by equation (13).
20
40
60
80
V. E XPERIMENTAL RESULTS
100
120
We experiment the proposed segmentation method using as filter response histograms, a set of co-occurrence distributions [27] computed for a displacement of one pixel in the four [2' H' 9[H'5 main directions ( ), a set of 4 distributions of the magnitude of Gabor filter responses, computed for several parameterizations (six radial frequencies I , and five
140
160 50
100
150
200
Fig. 3. Z EBRA IMAGE UNSUPERVISED SEGMENTATION .
[ 0'
2' 2'0[H' * !
r
orientations: and a set of distributions of the energy of the image wavelet packet coefficient computed for different bands for Haar wavelet. The computation of the co-occurrence [ distributions is issued from a quantization of the [9[ images into gray levels, where bins are exploited for Gabor and wavelet-based distributions. As most common methods, level-set functions are chosen to be the signed Euclidean distance to their zero level-sets. They are updated using gradient minimization techniques and re-initialized using a Partial Derivative Equations (PDE) based approach [21]. We initialize ` r according to an initial k-means segmentation (see section C). Figure and figure present the segmentation of a synthetic mosaic with three and respectively five textures selected from the Brodatz database.:The segmentation error rate are respectively and . In figure * we show the segmentation of a natural zebra image and in figure and we present the segmentations of real sidescan sonar images [28]. VI. C ONCLUSION We proposed a unsupervised criterion for unsupervised texture segmentation. The segmentation algorithm deals with multi-class case and it is based on texture characterization through marginal distributions of their filter responses. To improve computational efficiency and segmentation results, our
implementation includes an automatic feature selection step used to readjust the weights attached to each feature in the curve evolution equation that drives the segmentation. The main difference of this approach with texture unsupervised existing methods is the criterion formulation at region level. The segmentation algorithm is validated on several Brodatz texture collages and on real sonar images. R EFERENCES [1] L. Xiuwen and W. DeLiang, “Texture classification using spectral histograms,” IEEE Transactions on Image Processing, vol. 12, no. 6, pp. 661–670, 2003. [2] O.G Cula and K. Dana, “3d texture recognition using bidirectional feature histograms,” Proceedings of the 5th ACM SIGMM international workshop on Multimedia information retrieval, vol. 59, no. 1, pp. 33–60, 2003. [3] P. Nammalwar, O. Ghita, and P.F. Whelan, “Integration of feature distributions for color texture segmentation,” Conference on Pattern Recognition, vol. 1, pp. 716–719, 2005. [4] Q. Xu, J. Yang, and S. Ding, “Texture segmentation using lbp embedded region competition,” Electronic Letters on Computer Vision and Image Analysis, vol. 5, no. 1, pp. 41–47, 2004. [5] R. Fablet and P. Bouthemy, “Motion recognition using non parametric image motion models estimated from temporal and multiscale cooccurrence statistics,” Proceedings of the 5th ACM SIGMM international workshop on Multimedia information retrieval, vol. 59, no. 1, pp. 33–60, 2003. [6] S.C. Zhu, X.Liu, and Y.Wu, “Exploring texture ensembles by efficient markov chain monte carlo-toward a “trichromacy” theory of texture,” IEEE Transactions on Pattern Analysis and Machine Intelligence, vol. 22, no. 6, pp. 554–569, 2000.
50
100
150
200
250
50
100
150
200
250
300
Fig. 4. U NSUPERVISED SEGMENTATION OF 2 SEA - FLOOR TYPE SIDESCAN SONAR IMAGE (P ROJECT R EBENT, I FREMER ).
50
100
150
200
250
20
40
60
80
100
120
140
160
180
200
Fig. 5. U NSUPERVISED SEGMENTATION OF 3 SEA - FLOOR TYPE SIDESCAN SONAR IMAGE (P ROJECT R EBENT, I FREMER ).
[7] G.L. Gimel’farb, “Texture modeling by multiple pairwise pixel interactions,” IEEE Transactions on Pattern Analysis and Machine Intelligence, vol. 18, no. 11, pp. 1110–1114, 1996. [8] I. Karoui, R. Fablet, J.M. Boucher, and J.M. Augustin, “Regionbased image segmentation using texture statistics and level-set methods,” ICASSP, vol. 2, pp. 693–696, 2006. [9] G. Palubinskas, X. Descombes, and F. Kruggel, “An unsupervised clustering method using the entropy minimization,” International Conference on Pattern Recognition, vol. 2, pp. 1816–1818, 1998. [10] A. Lorette, X. Descombes, and J. Zerubia, “Fully unsupervised fuzzy clustering with entropy criterion,” International Conference on Pattern Recognition, vol. 3, pp. 986–989, 2000. [11] J. Kim, J.W. Fisher, A. Yezzi, M. Cetin, and A.S. Willsky, “A nonparametric statistical method for image segmentation using information theory and curve evolution,” IEEE Transaction on Image Processing, vol. 14, no. 10, pp. 1486–1502, 2005. [12] A. Herbulot, S. Jehan-Besson, M. Barlaud, and G. Aubert, “Shape gradient for multimodal image segmentation using mutual information,” International Conference on Image Processing, vol. 4, pp. 2729–2732, 2004. [13] M. Heiler and C. Schnorr, “Natural image statistics for natural image segmentation,” Proc. ICCV, vol. 2, pp. 1259–1266, 2003. [14] T. Brox, M. Rousson, R. Deriche, and J. Weickert, “Unsupervised segmentation incorporating colour, texture, and motion,” Computer Analysis of Images and Patterns. Lecture Notes in Computer Science, vol.
2756, pp. 353–360, 2003. [15] S. Kullback, Information theory and statistics, Wiley, New York, 1959. [16] T. Randen and J.H. Husoy, “Filtering for texture classification: a comparative study,” IEEE Transactions on Pattern Analysis and Machine Intelligence, vol. 21, no. 4, pp. 291–310, 1999. [17] E. Parzen, “On the estimation of a probability density function and the mode,” Annals of Mathematical Statistics, vol. 33, pp. 1065–1076, 1962. [18] I. Karoui, R. Fablet, J.M. Boucher, and J.M. Augustin, “Separabilitybased kullback divergence weighting and filter selection for texture classification and segmentation,” ICPR, pp. 18–21, 2006. [19] S. Osher and J. Sethian, “Fronts propagating with curvature dependent speed: Algorithms based on the hamiltonian-jacobi formulation,” Journal of computational physics, vol. 79, pp. 12–49, 1990. [20] I. Karoui, Segmentation par m´ethodes markoviennes et variationnelles des images textur´ees: application a` la caract´erisation sonar des fonds marins, PhD thesis, Ecole Nationale de T´el´ecommunocation de Bretagne, Brest, France. [21] C. Samson, L. Blanc-F´eraud, G. Aubert, and J. Zerubia, “A level set method for image classification,” International Journal on Computer Vision, vol. 40, no. 3, pp. 187–197, 2000. [22] H.K. Zhao, T. Chan, B. Merriman, and S. Osher, “A variational level set approach to multiphase motion,” J. Comp. Phy., vol. 127, pp. 179–195, 1996. [23] L. Vese and T. Chan, “A multiphase level set framework for image segmentation using the mumford and shah model,” International Journal of Computer Vision, vol. 50, no. 3, pp. 272–293, 2002. [24] I. Ben Ayed, A. Mitiche, and Z. Belhadj, “Multiregion level-set partitioning of synthetic aperture radar images,” IEEE Transactions On Pattern Analysis and Machine Intelligence, vol. 27, no. 5, pp. 793–800, 2005. [25] J.F Aujol, G. Aubert, and L. Blanc-F´eraud, “Wavelet-based level set evolution for classification of textured images,” IEEE Transactions on Image Processing, vol. 12, no. 12, pp. 1634–1641, 2003. [26] S. Jehan-Besson, M. Barlaud, and G. Aubert, “Image segmentation using active contours: calculus of variations for shape gradients?,” SIAM Journal APPL. MATH, vol. 63, no. 6, pp. 2128–2154, 2003. [27] R.M. Haralick, “Statistical and structural approaches to texture,” Proceedings of the IEEE, vol. 67, no. 5, pp. 786–804, 1979. [28] A. Ehrhold, D. Hamon, and B. Guillaumont, “The rebent monitoring network, a spatial integrated acoustic approach to survey nearshore macrobenthic habitats: application to the bay of concarneau (south brittany, france),” ICES Journal of Marine Science, vol. 63, pp. 1604– 1615, 2006.