Seabed segmentation using optimized statistics of sonar textures

IEEE TRANSACTIONS ON GEOSCIENCE AND REMOTE SENSING, VOL. 46, NO. 12, DECEMBER 2008 2 Seabed Segmentation Using Optimized Statistics of Sonar Textures 3 Imen Karoui, Ronan Fablet, Jean-Marc Boucher, Member, IEEE, and Jean-Marie Augustin 1 Abstract—In this paper, we propose and compare two supervised algorithms for the segmentation of textured sonar images, with respect to seafloor types. We characterize seafloors by a set of empirical distributions estimated on texture responses to a set of different filters. Moreover, we introduce a novel similarity measure between sonar textures in this feature space. Our similarity measure is defined as a weighted sum of Kullback–Leibler divergences between texture features. The weight setting is twofold. First, each filter is weighted according to its discrimination power: The computation of these weights are issued from a margin maximization criterion. Second, an additional weight, evaluated as an angular distance between the incidence angles of the compared texture samples, is considered to take into account sonar-image acquisition process that leads to a variability of the backscattered value and of the texture aspect with the incidence-angle range. A Bayesian framework is used in the first algorithm where the conditional likelihoods are expressed using the proposed similarity measure between local pixel statistics and the seafloor prototype statistics. The second method is based on a variational framework as the minimization of a region-based functional that involves the similarity between global-region texture-based statistics and the predefined prototypes. IE E Pr E oo f 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 1 Fig. 1. Typical sidscan sonar image (Rebent, IFREMER). 26 Index Terms—Active regions, angular backscattering, feature 27 selection, level sets, maximum marginal probability (MMP), seg28 mentation, sonar images, texture. I. I NTRODUCTION 29 COUSTIC remote sensing, such as high-resolution multibeam and sidescan sonars, provides new means for in situ observation of the seabed. The characterization of these high33 resolution sonar images is important for a number of practical 34 applications such as marine geology, commercial fishing, off35 shore oil prospecting, and drilling [1]–[4]. 36 The segmentation and the classification of sonar images with 37 respect to seafloor types (rocks, mud, sand, . . .) is the key 38 goal behind the analysis of these acoustic images. This task, 39 however, raises two major difficulties. The first task is to deal 40 with texture in these images. Previous methods are generally 41 based only on backscattered (BS) intensity models, and several 42 parametric families of probabilistic distribution functions have 30 31 32 A Manuscript received January 4, 2008; revised May 20, 2008 and August 11, 2008. I. Karoui, R. Fablet, and J.-M. Boucher are with the Department of Signal and Communications, École Nationale Supérieure des Télécommunications de Bretagne, 29238 Brest Cedex 3, France (e-mail: [email protected]). J.-M. Augustin is with the Department of Acoustics and Seismics, Institut Français de Recherche et d’Exploitation de la Mer, 29280 Brest, France. Color versions of one or more of the figures in this paper are available online at http://ieeexplore.ieee.org. Digital Object Identifier 10.1109/TGRS.2008.2006362 Fig. 2. BS evolution with the incidence angles for the three seafloor types of Fig. 1: sand, mud, and sand ripples. been suggested (Rayleigh distribution, K distribution, Weibull 43 distribution, etc.) [3]–[8]. These first-order statistics are not 44 sufficient when high-resolution sonar images involve textures, 45 which is the case of most sonar images (Fig. 1). 46 The other important issue arising in seabed texture charac- 47 terization is a built-in feature of sonar observation: The value 48 of BS measure depends both on the seafloor type and on the 49 incident angle of the reflected acoustic signal, ranging typically 50 from −85◦ to +85◦ . Fig. 2 shows the BS evolution for three 51 different seafloor types: sand, sand ripples, and mud. In addition 52 to the BS variability within incidence angles, seafloor textures 53 0196-2892/$25.00 © 2008 IEEE 2 IEEE TRANSACTIONS ON GEOSCIENCE AND REMOTE SENSING, VOL. 46, NO. 12, DECEMBER 2008 Fig. 3. Sonar image composed of sand ripples. (Left) [5◦ , 40◦ ]. (Right) [80◦ , 85◦ ]. 55 56 57 for two angular sectors. (Left) [5◦ , 40◦ ]. segmentation schemes. We first state the segmentation issue 93 as a Bayesian pixel-based labeling according to local texture 94 features. The second approach relies on a region-level vari- 95 ational framework, which resorts to a level-set minimization 96 of an energy criterion involving global-region-based seafloor 97 statistics. 98 This paper is organized as follows. The seafloor similarity 99 measure is introduced in Section II. The Bayesian segmentation 100 method is detailed in Section III. The region-based segmen- 101 tation algorithm is described in Section IV. Experiments are 102 reported and discussed in Section V, and conclusions are drawn 103 in Section VI. 104 IE E Pr E oo f are dependent on the incidence angles. Figs. 3 and 4 show a sonar image composed by sand ripples and rocks, respectively, for two angular sectors: [80◦ , 85◦ ] and [5◦ , 40◦ ]. The texture of sand ripples shows a loss of contrast in the specular sec58 tor [5◦ , 40◦ ]: The steep grazing angle reduces the backscatter 59 differences between facing and trailing slopes, while at low 60 incidence angles, much of the variation is lost due to increasing 61 sonar shadow. A similar loss of contrast is observable in the 62 rock samples (Fig. 4). BS behavior according to the incidence 63 angles has been of wide interest for sonar imaging [3], [9]–[12]. 64 Parametric and nonparametric techniques have been proposed 65 to model sonar-image behavior with respect to the incidence 66 angle variations. The effect of the incidence angle on the BS 67 has also been explored as a discriminating feature for seafloor 68 recognition [4], [13], [14]. However, no studies have proposed 69 a method to accurately compensate this phenomena because 70 of the joint dependence of the seafloor types and on the local 71 bathymetry which is generally unknown for sidescan sonar 72 images. To our knowledge, the effects of the incidence angle on 73 textured seabed features have not been addressed for segmen74 tation issues. Only some studies were interested in simulating 75 the behavior of oriented and textured seafloor types [15]–[18]. 76 These methods are mainly based on shape-from-shading [19] 77 and were restricted to synthetic images or to real sonar images 78 involving only one seafloor type. 79 In this paper, we aim at using texture information within 80 sonar seabed images and at developing segmentation algo81 rithms that take into account angular variability of BS and 82 textural features. We propose to characterize seafloors by a 83 wide set of marginal distributions of their filter responses, and 84 we measure seafloor similarities according to a weighted sum 85 of Kullback–Leibler divergences [28] in this feature space. 86 To cope with seafloor angular dependence, we introduce an 87 additional weighting factor, evaluated as an angular distance 88 between the compared texture samples: This angular distance 89 is measured according to Gaussian kernels, whose variance 90 sets the level of the angular variability, depending on textures 91 and seafloor types. The proposed incidence-angle-and-texture92 based similarity measure is exploited to develop two different 54 Fig. 4. Rock texture (Right) [80◦ , 85◦ ]. II. S ONAR -T EXTURE S IMILARITY M EASURE 105 Texture-based segmentation of seabed sonar images gener- 106 ally relies on Haralick [30] parameters or scalar spectral and 107 filter coefficients to model textures [20]–[22]. Recently, in the 108 field of texture analysis, features computed as statistics of local 109 filter responses have been shown to be relevant and discriminant 110 texture descriptors [23]–[27]. Motivated by these studies, we 111 propose to use texture features computed as marginal distrib- 112 utions of a wide set of filter seafloor responses. Each seafloor 113 type denoted by Tk is characterized by a set Qk composed of 114 the marginal distributions of the seafloor with respect to the 115 predefined filters. The following conditions are issued. 116 1) 121 cooccurrence distributions with horizontal and ver- 117 tical displacements denoted by dx and dy , respectively, 118 (dx , dy ) ∈ {0, 1, . . . , 10}. 119 2) 50 distributions of the magnitude of Gabor filter 120 responses, computed for combination of parameters 121 (f0 , σ, θ), where f0 is the radial frequency, σ is 122 the standard√deviation, and θ is the √ orientation, such 123 that f0 ∈ { 2/2k }k=1:6 , σ ∈ {2k 2}k=2:5 , and θ ∈ 124 {0◦ , 25◦ , 45◦ , 90◦ , 135◦ }. 125 3) 48 distributions of the energy of the image wavelet packet 126 coefficient computed for different bands (we used three 127 wavelet types: Haar, Daubechies, and Coiflet). 128 KAROUI et al.: SEABED SEGMENTATION USING OPTIMIZED STATISTICS OF SONAR TEXTURES 3 Gaussian mixture, and the assignment likelihood p(Θj /zs ) is 147 given by 148   −(Θj −θs )2 exp σ2 .  j p(Θj /zs ) =  (1) −(Θi −θs )2 J exp 2 i=1 σ i Hence, the assignment of a given region R to an angular 149 subdomain is given by 150  1 p(Θj /zs ), ds. (2) πR (j) = |R| s IE E Pr E oo f Fig. 5. Cooccurrence distributions computed for parameters (d, θ) = (6, 0◦ ). (Left) [5◦ , 40◦ ]. (Right) [80◦ , 85◦ ]. The angular bounds, defined by the parameters Θj and σj , 151 are set experimentally, and each seafloor type is character- 152 ized by Qk = {Qf,j (Tk )}f =1:F,j=1:J , its filter responses es- 153 timated on the J angular domains. Note that f accounts 154 both for filter types and various associated parameteri- 155 zations: f = 1 → 121 refers to cooccurrence distributions, 156 f = 122 → 171 to Gabor filter, and f = 172 → 219 to 157 wavelet-based distributions. These statistics are computed us- 158 ing Parzen window estimation [29]. Formally, Qf,j (Tk , α) = 159 (1/πTk (j)) Tk p(Θj /zs )gσf (hf (s) − α)ds, where hf is the 160 filter response indexed by f (for a cooccurrence matrix with pa- 161 rameters μ = (d, θ), hf : Ω → [1, N g] × [1, N g], and hf (s) = 162 (I(s), I(s + μ)), where I(s) is the gray value at pixel s and N g 163 is the total gray level number) and gσf is a Gaussian kernel with 164 zero mean and standard deviation σf . 165 We define the similarity measure between a texture sample 166 T with feature set D(T ) = {Df,j (T )}f =1:F,j=1:J and a given 167 seafloor type Tk as follows: 168 J  k  KLΘ w Q , D(T ) = Fig. 6. Cooccurrence distributions computed for parameters (d, θ) = (Left) [5◦ , 40◦ ]. (Right) [80◦ , 85◦ ]. (6◦ , 0). Sonar-texture variability with respect to the incidence angles induces a variability in texture features. Figs. 5 and 6 show the cooccurrence matrices of the images involving sand ripples 132 and rocks shown in Figs. 3 and 4, respectively. A change in 133 cooccurrence distributions can be noticed. In the angular sector 134 [5◦ , 40◦ ], the cooccurrence is bimodal due to the alternation of 135 dark and light values in the image of sand ripples, whereas in 136 the sector [80◦ , 85◦ ], the loss of contrast between dark and light 137 values in the image of sand ripples leads to unimodal distribu138 tion. For cooccurrence distributions related to rock samples, a 139 change in the variance is depicted. 140 To deal with this problem, we propose to define angular 141 subdomains, in which texture characteristics can be regarded 142 as homogeneous. Each angular sector is indexed by j and 143 is characterized by a mean incidence-angle value Θj . For144 mally, this leads in introducing a state variable zs assign145 ing pixel s to a given angular subdomain. (θs , zs ), where 146 θs is the incidence angle of pixel s, is then modeled as a 129 130 131 F j=1 f =1   πT (j)wf2 KL Qkf,j , Df,j (T ) (3) where {wf2 }f =1:F are the feature weights such that 169 F 2 f =1 wf = 1, Df,j (T ) is the distribution related to 170 filter f and estimated on  region T for angular sector j: 171 Df,j (T, α) = (1/πT (j)) T p(Θj /zs )gσf (hf (s) − α)ds, and 172 KL denotes for the Kullback–Leibler divergence [28]; 173 for  two probability distributions Q and D, KL(Q, D) = 174 Q(α)log(Q(α)/D(α))dα. 175 The resulting weighting factors {wf2 }f =1:F are exploited on 176 the one hand for filter selection to keep only the distributions 177 corresponding to the significant weights and, on the other hand, 178 for the definition of an optimized texture-based similarity mea- 179 sure KLΘ w given the selected distributions. In a supervised con- 180 text, the weights are estimated from a training set T composed 181 of N -labeled texture samples: T = {(D(T ), sT )}, where sT is 182 the class of sample T . Formally, {wf2 }f =1:F are issued from the 183 maximization of the global margin defined as follows: 184 MwT = Mw (T ) (4) T ∈T where  d  Θ sT T Mw (T ) = KLΘ w Q , D(T ) − KLw (Q , D(T )) 185 (5) 4 186 187 IEEE TRANSACTIONS ON GEOSCIENCE AND REMOTE SENSING, VOL. 46, NO. 12, DECEMBER 2008 where dT is the nearest class to T different from sT according to the similarity measure KLΘ w  k  (6) dT = arg min KLΘ w Q , D(T ) . k=sT iterative conditional estimation procedure [42]. Using Bayes 227 rule, the posterior distribution is expressed as follows: 228 p(X, Y ) = p(Y /X)p(X) = p(X) The maximization of the margin criterion MwT is carried out using a stochastic gradient method as detailed in our previous 190 work [39]. 188 189 191 III. B AYESIAN S ONAR -I MAGE S EGMENTATION As far as Bayesian image segmentation is concerned, the most popular criteria are the maximum a posteriori (MAP) 194 [31] and the maximum marginal probability (MMP) [32]. It 195 has been shown that the MMP estimation criterion is more 196 appropriate for image segmentation than the MAP criterion 197 [32]. The MAP estimate assigns the same cost to every incorrect 198 segmentation regardless the number of pixels at which the 199 estimated segmentation differs from the true one, whereas the 200 MMP algorithm minimizes the expected value of the number of 201 misclassified pixels. As shown in [32], the MMP procedure is 202 equivalent to maximizing the marginal of the class labels. Let 203 us introduce the following notations: 204 1) S, the image lattice composed of N pixels; 205 2) X = {xs }s∈S , the label random field; 206 3) Y = {ys }s∈S , the random field of the observations, the 207 textural feature in our case. 208 Formally, the MMP scheme is equivalent to the maximi209 zation of x̂MMP = arg max p(xs = k/y). s k∈{1:K} (7) In general, pixel conditional likelihood is computed according to local texture features like Gabor or wavelet coefficients. Here, we aim at using the proposed similarity measure KLΘ w. We associate to each pixel s a set of features D(Ws ) = 214 {Df (Ws )}f =1:F estimated according to a Parzen estimation 215 method [29], within a square window Ws centered at s. The 216 window size that we denote by TW is set by the user according 217 to texture coarseness. In our case, the observation denoted 218 by a random field Y is specified by {ys = D(Ws )}s∈S . The 219 conditional likelihood at each pixel s with respect to class k is 220 then defined from the similarity measure KLΘ w by 210 211 212 213 k −KLΘ w (Q ,D(Ws )) 221 222 223 exp . p(ys /xs = k) = K i −KLΘ w (Q ,D(Ws )) i=1 exp s∈S t∈cs −αc (1−δ(xs ,xt ))+log(p(ys /xs )) . (10) The maximization of local probabilities p(xs = k/y) is car- 229 ried out as follows [32]: We use the Gibbs sampler to gen- 230 erate a discrete-time Markov chain X(t) which converges in 231 distribution to a random field with probability mass functions 232 P (X/Y ). The marginal conditional distributions p(xs = k/y) 233 are then approximated as the fraction of time the Markov chain 234 spends in state k for each class k and each pixel s. The MMP 235 segmentation steps are the following. 236 1) Simulation of Tmax realizations of x1 , x2 , . . . , xTmax of 237 X using Gibbs sampler. 238 2) Using the realizations x1 , x2 , . . . , xTmax , p(xs /y) is esti- 239 mated using the frequency of each realizations 240     δ x1s − k + · · · + δ x1s − k p(xs = k/y) = . Tmax s∈S t∈cs αc (1 − δ(xs , xt )) 3) Choose as xs the class that maximizes p(xs = k/y). 241 IV. V ARIATIONAL S ONAR -I MAGE S EGMENTATION 242 Unlike the Bayesian scheme, the second approach is stated 243 at a region level as the minimization of a constrained energy 244 criterion E({Ωk }k=1:K ) = E1 + E2 , where Ωk is the domain 245 composed of all pixels attributed to the class k, E1 is a texture- 246 based data-driven term, and E2 is a regularization term detailed 247 as follows. 248 A. Functional Terms E ({Ωk }k=1:K ) is evaluated as the log-likelihood of a given 250 partition with respect to texture models. It is evaluated as the 251 sum of the similarities according to the measure KLΘ w between 252 each region Ωk and its corresponding class Tk 253 K E1 ({Ωk }k=1:K ) = k=1  k  KLΘ w Q , D(Ωk ) (11) where Df,j (Ωk ) is the marginal distribution of the image 254 response to the filter indexed by f . For angular domain j, 255 Df,j (Ωk ) is estimated according to Parzen method [29]. 256 E2 is a regularization term, it penalizes the lengths of region 257 contours and is expressed by 258 (9) where δ is the delta function and αc ∈ {αH , αV , αD } are real parameters assigned, respectively, to horizontal, vertical, and 226 diagonal cliques. These parameters are estimated using the 249 1 (8) As a prior PX , we consider a Markov random field associated to an eight-neighborhood system with potential functions given by U2 (x) = 224 225 ≈ exp   IE E Pr E oo f 192 193 p(ys /xs ) s K E2 = k=1 γk |Γk |, γk ∈ ℜ+ (12) where |Γk | is the length of the contour Γk associated to the 259 region Ωk . 260 KAROUI et al.: SEABED SEGMENTATION USING OPTIMIZED STATISTICS OF SONAR TEXTURES 261 B. Computation of the Evolution Equation We solve for the minimization of the functional E using a gradient-descent technique. It relies on the computation of 264 the first derivative of E according to regions {Ωk }k=1:K . The 265 evolution equation of region contours {Γk }k=1:K is then given 266 by the following dynamic scheme [33]: 262 263 ∂Γk (x,t) = Fk (x, t)Nk ∂t Γk (x, 0) = Γ0k (13) 5 where ∂ϕik /∂t, i = 1, 2, 3 are the evolution equation terms associated, respectively, to functionals Ei , i = 1, 2, 3. The derivatives of the energy terms E2 and E3 are directly estimated from level-set functions [35], [36] 287 288 289 290 ∇ϕk ∂ϕ2k (s, t) = γk δα (ϕk )div ∀k ∈ {1 : K} (19) ∂t |∇ϕk | K  ∂ϕ3k (s, t) =−δα (ϕk )λ (Hα (ϕk )−1) ∀k ∈ {1 : K}. ∂t k=1 where Nk is the unit inward normal to Γk at pixel x and at time t and Fk is the velocity field (in our case, Fk = ∇Ek , the 269 derivative of E with respect to Γk ). 270 The explicit implementation of the curve evolution according 271 to the latter dynamic scheme using a difference-approximation 272 scheme cannot deal with topological changes of the moving 273 front. This could be avoided by introducing the level-set method 274 proposed by Osher and Sethian [34]. The basic idea of the 275 method is the implicit representation of the moving interface Γ 276 by a higher dimensional hypersurface ϕ (the level-set function) 277 such that the zero level set of ϕ is actually the set of Γ and 267 268 K γk lim k=1 281 ∀s ∈ Ωk . (14) E2 can be expressed using level-set functions ϕk [36] E2 = α→0  δα (ϕk )|∇ϕk |ds (15) Ω where δα is a regularized delta function    1 1 + cos πs , 2α α δα (s) = 0, if |s| ≤ α if |s| < α. (16) In order to cope with multiclass segmentation and to fulfill the image-partition constraint, we use an additional term E3 given 284 by the following functional [38]: 282 283 λ E3 = 2   Ω 285 286 K k=1 F f =1    k wf2 p(Θj /zs ) KL Qkf,j , Df,j  − Qkf,j /Df,j ∗ gσf (hf (s)) − 1 |∇ϕk | The evolution of the contours {Γk }k=1:K (13) is then equivalent to the evolution of the level-set functions {ϕk }k=1:K [34] ∂ϕk (s, t) = Fk (s, t) |∇ϕk (s, t)| ∂t 280 J ∂ϕ1k (s, t) = − ∂t j=1 IE E Pr E oo f 278 The evolution equation related to E1 is more complex, since it 291 involves computations over the spatial support of each region. 292 To differentiate E1 , we use shape-derivative tools, particularly 293 the Gâteaux derivative theorem given in [37]. As detailed in the 294 Appendix, it leads to 295  Ωinside = {s ∈ Ω/ϕ(s) > 0} Ωoutside = {s ∈ Ω/ϕ(s) < 0} . 279 (20) 2 Hα (ϕk ) − 1 ds, λ ∈ ℜ+ where Hα is a regularized Heaviside function ⎧    ⎨ 12 1 + αs + π1 sin πs , if |s| ≤ α α Hα (s) = 1, if s > α ⎩ 0, if s < −α. As E = E1 + E2 + E3 , we have ∂ϕk ∂ϕ1k ∂ϕ2k ∂ϕ3k (s, t) = (s, t) + (s, t) + (s, t) ∂t ∂t ∂t ∂t (17) (21) where ∗ is the convolution symbol. 296 The evolution equation related to the energy E1 (21) is 297 composed of two terms. 298 J  F 2 k 1) A global term − j=1 f =1 wf p(Θj /zs )KL(Qf,j , 299 Df,j ): This term is always negative or null. It is a contrac- 300 tion force that  reducesthe size of heterogeneous regions. 301 F J 2 k 2) A local term f =1 wf p(Θj /zs )((Qf,j /Df,j ) − 302 j=1 1) ∗ gσf (hf (s)): This term locally compares the features 303 values at each pixel. This term can be positive or neg- 304 ative and aims at readjusting the statistics inside the 305 regions Df,j (Ωk ) to fit to prototype models Qkf,j . The 306 contribution of each descriptor f is weighted by wf2 and 307 p(Θj /zs ), the relative contribution of descriptor f and of 308 angular sector j. 309 The overall evolution equations of the contours {Γk }k=1:K 310 are the following: 311 ∂ϕk (s, t) ∂t = −δα (ϕk ) ⎡ J ×⎣ F wf2 p(Θj /zs ) j=1 f =1 ×  KL(Qkf,j , Df,j )− (18) + γk div ∇ϕk −λ |∇ϕk |  Qkf,j ∗ gσf (hf (s))+1 Df,j  K (Hα (ϕk )−1) k=1 ∀k. (22) We apply these coupled evolution equations until convergence. 312 6 IEEE TRANSACTIONS ON GEOSCIENCE AND REMOTE SENSING, VOL. 46, NO. 12, DECEMBER 2008 Fig. 8. Feature weights. 1 → 121 cooccurrence distributions, 122 → 171 Gabor filter-based features and 172 → 219 wavelet-based distributions. (a) I1. (b) I5. IE E Pr E oo f TABLE I SEGMENTATION ERROR RATES Fig. 7. Test images and their manual segmentation (in black line). 313 V. E XPERIMENTAL R ESULTS AND D ISCUSSION In previous work, we have tested the method on various optic textures (Brodatz textures). The method was compared 316 to other texture-classification methods, and some results are 317 reported in [39] and [40]. Here, we evaluate the proposed 318 seabed-segmentation technique for different real sonar images 319 acquired by a sidescan sonar, as part of a natural seabed320 mapping project (IFREMER, Rebent Project) [41]. A reference 321 interpretation by an expert is available [41]. Fig. 7 shows the 322 set of images on which we carried out the experiments. We 323 superimposed on these images the manual expert segmentation. 324 Image I1 is composed of three seafloor types: rock, mud, and 325 marl ripples [41]; I2 of rock, marl ripples, and mud seafloors; 326 I3 of mud, sand, and marl ripples; I4 of marl and marl ripples; 327 and I5 of sand, sand ripples, and rock. For I2 and I3, the angular 328 variability of the sefloors, particularly marl ripples (for I2) and 329 marl (for I3), is visually clear. 330 For all these images, we first determine the most discriminant 331 features among the initial set of 219 features: We apply the 332 algorithm described in Section II and detailed in [39], and 333 we keep only the feature set such that the cumulative sum 334 of weights exceeds 0.9. Only a small number of features are 335 retained. For example, Fig. 8 shows the plot of feature weights 336 computed for image I1 and I5. For I1, the two cooccurrence 337 matrices account for more than 90% of the total weight sum; 338 these cooccurrence are computed for parameters (dx , dy ) ∈ 339 {(1, 4), (2, 1)}. For I5, the cooccurrence distributions com340 puted for parameters (dx , dy ) ∈ {(2, 1), (2, 2)} are selected. 341 For sonar images, we noticed that cooccurrence matrices are the 314 315 most selected features. In previous work on Brodatz textures 342 [39], [40], we remarked that Gabor and wavelet filters were 343 selected for oriented textures, whereas cooccurrence distribu- 344 tions, which, in addition to the detection of texture structures, 345 detect the intensity values change, are selected in the case of 346 texture having different intensity values and texture with regular 347 motifs. 348 For the five test images, three segmentation algorithms are 349 compared. 350 1) The maximum-likelihood (ML) segmentation denoted by 351 ML. This method consists in maximizing at each pixel the 352 conditional probability p(ys /xs ) given by (8) 353 xˆs = arg max p(ys /xs = k). k (23) Several analysis-window sizes are compared: TW ∈ {7 × 354 7, 17 × 17, 33 × 33}. 355 2) The MMP segmentation described in Section III, applied 356 for several analysis-window sizes: TW ∈ {7 × 7, 17 × 357 17, 33 × 33}. 358 3) The region-based variational segmentation described in 359 Section IV that we denote by V ar. 360 Table I summarizes the different error classification rates for all 361 segmentations. 362 All segmentation methods give quite good results according 363 to the mean classification error rates. MMP and variational 364 approaches are more efficient than the ML-based segmentation 365 because they take into account the spatial dependence between 366 pixels. The difference between these later approaches (MMP 367 and variational methods) mainly lies in the accuracy of the 368 localization of region boundaries and in the dependence of 369 MMP-based segmentation on the window size TW . In fact, 370 for the MMP segmentation, the use of small neighborhood 371 (TW = 7 × 7) leads to more accurate region frontiers but 372 small misclassified patches appear because of the neighborhood 373 KAROUI et al.: SEABED SEGMENTATION USING OPTIMIZED STATISTICS OF SONAR TEXTURES MMP segmentation of I1 using TW = 7 × 7, τ = 12.5%. Fig. 10. MMP segmentation of I1 using TW = 17 × 17, τ = 9.2%. inefficiency for texture characterization. Conversely, the use of large window sizes (TW = 33 × 33) resorts to a lack of 376 accuracy in the localization of the boundaries of the seabed 377 regions because texture features extracted for pixels close to 378 region boundaries involve a mixture of texture characteris379 tics. The variational region-based approach does not need the 380 choice of an analysis window and operates globally on region 381 composed of pixels belonging to the same class. It resorts 382 to a tradeoff between segmentation accuracy and region ho383 mogeneity. Figs. 9–13 show examples of the dependence of 384 MMP segmentation on the sizes of analysis window and on the 385 robustness of the variational approach. In Fig. 14, we plot the 386 mean error rate for different window sizes. 387 We note that MMP and variational segmentations give sim388 ilar results when the analysis-window size is well chosen for 389 instance, TW = 17 × 17 for image I2 (Figs. 15 and 16), but 390 its performance depends a lot on the choice of this parameter. 391 This method can however be appropriate when the aim of the 392 segmentation is to detect texture regions and without seeking 393 for accurate boundaries. 394 The variational approach is also interesting because it is 395 much faster than the MMP segmentation. Being deterministic, 396 the variational approach can be very fast particularly if we use 397 appropriate initialization such as an initial segmentation based 398 on the ML criterion. Whereas the MMP segmentation needs 374 375 Fig. 11. MMP segmentation of I1 using TW = 33 × 33, τ = 13%. IE E Pr E oo f Fig. 9. 7 Fig. 12. Region-based segmentation of I1, τ = 7.5%. a large number of iterations, each iteration is also complex 399 and involves Gibbs sampling. For our implementations, the 400 convergence time of a variational image typically corresponds 401 to one iteration of the MMP algorithm. 402 To stress the interest of taking into account texture variability 403 with respect to incidence angles, additional segmentation re- 404 sults are reported for image I2 and I3 for which the seafloor 405 texture variability is clearer. In Table II, we summarize the 406 segmentation error rates for the segmentation with (J = 3) and 407 with no angular weighting, i.e., using only one angular domain 408 J = 1, for ML, MMP, and variational methods. In Figs. 17 and 409 18 are reported the associated segmentation results. It can be 410 noticed that a classical segmentation (without taking into ac- 411 count texture variability within incidence angles) cannot distin- 412 guish between visually similar seafloors (mud and marl ripples 413 for I2 and marl and sand for I3) near the specular domain. 414 VI. C ONCLUSION 415 We proposed two segmentation algorithms for sonar-image 416 segmentation: a Bayesian algorithm using local statistics and a 417 region-based variational algorithm, both based on a novel sim- 418 ilarity measure between seafloor-type images according to the 419 statistics of their responses to a large set of filters. This similar- 420 ity measure is expressed as a weighted sum of Kullback–Leibler 421 8 IEEE TRANSACTIONS ON GEOSCIENCE AND REMOTE SENSING, VOL. 46, NO. 12, DECEMBER 2008 IE E Pr E oo f Fig. 13. Segmentations of I4. (a) MMP: TW = 7, τ = 4%. (b) MMP: TW = 33, τ = 6%. (c) Variational method: τ = 3%. Fig. 16. Variational-based segmentation of I2, τ = 3%. Fig. 14. Mean segmentation error rate for several window sizes. TABLE II SEGMENTATION ERROR RATES Fig. 15. MMP-based segmentation of I2 TW = 17, τ = 4%. divergences between individual seafloor filter-response statistics. The resulting weighting factors are exploited on the one hand for filter selection and, on the other hand, for taking into account the incidence angular dependence of seafloor tex426 tures. The conclusion is the cooccurrence matrices outperform 427 the other features for our sonar images. The results show 428 that the performance of the Bayesian approach depends on 429 the size of the analysis window. For pixel-based segmenta430 tion (ML and MMP segmentations), the size must be tuned 431 according to the coarseness of given textures. The results 432 also stress the suitability of the region-based approach as 433 compared to the Bayesian pixel-based scheme for texture seg434 mentation and demonstrate the interest of taking into account 422 423 424 425 Fig. 17. I2 region-based segmentation with and without angular weighting. (a) J = 1, τ = 23.5%. (b) J = 3, τ = 3%. Fig. 18. I3 MMP-based segmentation with and without angular weighting. (a) J = 1, τ = 16%. (b) J = 3, τ = 8%. the angular backscatter variabilities to discriminate between 435 seafloor types particularly near the specular sector. 436 KAROUI et al.: SEABED SEGMENTATION USING OPTIMIZED STATISTICS OF SONAR TEXTURES A PPENDIX E VOLUTION E QUATION C OMPUTATION 437 438 441 Using the shape-derivative tools, we want to differentiate the functional  k  F (Ωk ) = KLΘ w Q , D(Ωk ) which can be written as follows: J F πΩk (j)wf2 j=1 f =1 442 443 444  Qkf,j (α) log Rf Qkf,j (α) Df,j (Ωk , α)  Then, we have J F wf2 πΩk (j) j=1 f =1 J F Finally, we get 457 ) dHjk (Ωk , V   = p(Θj /zs )  Γk Qkf,j (α) Qf (α) gσf (hf (s) − α) − G1 (f, α, Ωk ) πΩk (j) )  Hjk . Γk + πΩk (j) Hjk dα  )dα. dHjk (Ωk , V  · Nk )da(s) × (V   )dα dHjk (Ωk , V Rf 1 = πΩk (j) Rf    · Nk )da(s) × (V    Qkf,j (α) 1 = gσ (hf (s)−α) Qf (α) p(Θj /zs ) πΩk (j) Df,j (Ωk , α) f Rf  wf2 dπΩk (j)(Ωk , V j=1 f =1 Qkf,j ∂Hjk = . ∂πΩk (j) πΩk (j) Γk The Gâteaux derivative of F (Ωk ) in the direction of a vector  is given by field V )= dF (Ωk , V Qkf,j ∂Hjk =− ∂G1 G1 According to the previous theorem, we get 456   · Nk )da(s).  ) = − p(Θj /zs )gσ (hf (s) − α) (V dG1 (Ωk , V f p(Θj /zs )gσf (hf (s) − α) ds. F (Ωk ) = 446 dα. 2) Ωk 447  Let us introduce the following notations. 1)   Qkf,j (α) k k Hj = Qf,j (α) log . Df,j (Ωk , α) G1 = 445  455 Using the chain rule, we get   ∂Hjk ∂Hjk  )=  )+ ) dHjk (Ωk , V dG1(Ωk , V dπΩk(j)(Ωk , V ∂G1 ∂πΩk (j) IE E Pr E oo f 439 440 9  Γk   Qkf,j ∗ gσf (hf (s)) − 1 p(Θj /zs ) Df,j (Ωk )  · Nk )da(s). × (V (24) Finally, according to (24) 458 Rf 448 Theorem: See [37]. 449  The Gâteaux derivative 450 Ω k(s, Ω)ds is given by )= dK(Ω, V  ) dF (Ωk , V of a functional of the type K(Ω) =  )ds − ksh (s, Ω, V  N  )da(s) k(s, Ω)(V Γ Ω J =− F p(Θj /zs )wf2 j=1 f =1       k  Qkf,j × KL Qf,j , Df,j − −1 ∗ gσf (hf (s)) Df,j (Ωk ) Γk  ) is where Ω denoted a region and Γ its boundary; ksh (s, Ω, V the shape derivative of k(s, Ω). 453 In our case   ) = − p(Θj /zs )(V  · Nk )da(s) dπΩk (j)(Ωk , V 451 452 Γk 454 Hjk can be written as Hjk = Qkf,j (α) log  πΩk (j)Qkf,j (α) G1 (Ωk , f, j, α)  .    · Nk da(s). × V R EFERENCES 459 [1] P. Brehmer, F. Gerlotto, J. Guillard, F. Sanguinède, Y. Guénnegan, and D. Buestel, “New applications of hydroacoustic methods for monitoring shallow water aquatic ecosystems: The case of mussel culture grounds,” Aquat. Living Resour., vol. 16, no. 3, pp. 333–338, Jul. 2003. [2] M. Mignotte, C. Collet, P. Perez, and P. Bouthemy, “Hybrid genetic optimization and statistical model based approach for the classification of shadow shapes in sonar imagery,” IEEE Trans. Pattern Anal. Mach. Intell., vol. 22, no. 2, pp. 129–141, Feb. 2000. 460 461 462 463 464 465 466 467 10 [3] L. Hellequin, J. M. Boucher, and X. Lurton, “Processing of highfrequency multibeam echo sounder data for seafloor characterization,” IEEE J. Ocean. Eng., vol. 28, no. 1, pp. 78–89, Jan. 2003. [4] G. Le Chenadec and J. M. Boucher, “Sonar image segmentation using the angular dependence of backscattering distributions,” in Proc. Oceans—Europe, Brest, France, 2005, vol. 1, pp. 147–152. [5] G. Le Chenadec, J.-M. Boucher, and X. Lurton, “Angular dependence of K-distributed sonar data,” IEEE Trans. Geosci. Remote Sens., vol. 45, no. 5, pp. 1224–1236, May 2007. [6] E. Jakeman, “Non-Gaussian models for the statistics of scattered waves,” Adv. Phys., vol. 37, no. 5, pp. 471–529, Oct. 1988. [7] A. P. Lyons and D. A. Abraham, “Statistical characterization of highfrequency shallow-water seafloor backscatter,” J. Acoust. Soc. Amer., vol. 106, no. 3, pp. 1307–1315, Sep. 1999. [8] D. A. Abraham and A. P. Lyons, “Novel physical interpretations of k-distributed reverberation,” IEEE J. Ocean. Eng., vol. 27, no. 4, pp. 800– 813, Oct. 2002. [9] N. C. Mitchell, “Processing and analysis of Simrad multibeam sonar data,” Mar. Geophys. Res., vol. 18, no. 6, pp. 729–739, Dec. 1996. [10] G. Le Chenadec, J. M. Boucher, X. Lurton, and J. M. Augustin, “Angular dependence of statistical distributions for backscattered signals: Modeling and application to multibeam echosounder data,” in Proc. Oceans, Sep. 2003, vol. 2, pp. 897–903. [11] D. R. Jackson, D. P. Winbrenner, and A. Ishimaru, “Application of the composite roughness model to high-frequency bottom backscattering,” J. Acoust. Soc. Amer., vol. 79, no. 5, pp. 1410–1422, May 1986. [12] C. de Moustier and D. Alexandrou, “Angular dependence of 12-kHz seafloor acoustic backscatter,” J. Acoust. Soc. Amer., vol. 90, no. 1, pp. 522–531, Jul. 1991. [13] J. E. Hughes-Clarke, “Toward remote seafloor classification using the angular response of acoustic backscatter: A case study from multiple overlapping GLORIA data,” IEEE J. Ocean. Eng., vol. 19, no. 1, pp. 112– 127, Jan. 1994. [14] E. Pouliquen and X. Lurton, “Automated sea-bed classification system for echo-sounders,” in Proc. Oceans, Oct. 1992, vol. 1, pp. 317–321. [15] D. Langer and M. Herbert, “Building qualitative elevation maps from side scan sonar data for autonomous underwater navigation,” in Proc. IEEE Int. Conf. Robot. Autom., Apr. 1991, vol. 3, pp. 2478–2483. [16] R. Li and S. Pai, “Improvement of bathymetric data bases by shape from shading technique using side-scan sonar images,” in Proc. OCEANS, Oct. 1991, vol. 1, pp. 320–324. [17] E. Dura, J. Bell, and D. Lane, “Reconstruction of textured seafloors from side-scan sonar images,” Proc. Inst. Elect. Eng.—Radar Sonar Navig., vol. 151, no. 2, pp. 114–125, Apr. 2004. [18] A. E. Johnson and M. Hebert, “Seafloor map generation for autonomous underwater vehicle navigation,” Auton. Robots, vol. 3, no. 2, pp. 145–168, Jun. 1996. [19] R. Zhang, P. Tsai, J. E. Cryer, and M. Shah, “Shape from shading: A survey,” IEEE Trans. Pattern Anal. Mach. Intell., vol. 21, no. 8, pp. 690– 705, Aug. 1999. [20] T. H. Eggen, “Acoustic sediment classification experiment by means of multibeam echosounders,” Acoust. Classification Mapping Seabed, vol. 15, no. 2, pp. 437–444, 1993. [21] L. M. Linnett, S. J. Clarke, and D. R. Carmichael, “The analysis of sidescan sonar images for seabed types and objects,” in Proc. 2nd Eur. Conf. Underwater Acoust., 1994, vol. 22, pp. 187–194. [22] M. Lianantonakis and Y. R. Petillot, “Sidescan sonar segmentation using active contours and level set method,” in Proc. OCEANS—Europe, Brest, France, 2005, vol. 1, pp. 719–724. [23] L. Xiuwen and W. DeLiang, “Texture classification using spectral histograms,” IEEE Trans. Image Process., vol. 12, no. 6, pp. 661–670, Jun. 2003. [24] O. G. Cula and K. Dana, “3D texture recognition using bidirectional feature histograms,” Int. J. Comput. Vis., vol. 59, no. 1, pp. 33–60, Aug. 2004. [25] R. Fablet and P. Bouthemy, “Motion recognition using nonparametric image motion models estimated from temporal and multiscale co-occurrence statistics,” IEEE Trans. Pattern Anal. Mach. Intell., vol. 25, no. 12, pp. 1619–1624, Dec. 2003. [26] P. Nammalwar, O. Ghita, and P. F. Whelan, “Integration of feature distributions for colour texture segmentation,” in Proc. Int. Conf. Pattern Recog., Aug. 2005, vol. 1, pp. 716–719. [27] Q. Xu, J. Yang, and S. Ding, “Texture Segmentation using LBP embedded region competition,” Electron. Lett. Comput. Vis. Image Anal., vol. 5, no. 1, pp. 41–47, 2004. [28] S. Kullback, “On information and sufficiency,” in The Annals of Mathematical Statistics, vol. 22. New York: Wiley, 1951. [29] E. Parzen, “On estimation of a probability density function and mode,” Ann. Math. Stat., vol. 33, no. 3, pp. 1065–1076, Sep. 1962. [30] R. M. Haralick, “Statistical and structural approaches to texture,” Proc. IEEE, vol. 67, no. 5, pp. 786–804, May 1979. [31] S. Geman and G. Geman, “Stochastic relaxation, Gibbs distributions, and the Bayesian restoration of images,” IEEE Trans. Pattern Anal. Mach. Intell., vol. PAMI-6, no. 6, pp. 721–741, Nov. 1984. [32] J. Marroquin, S. Mitter, and T. Poggio, “Probabilistic solution of ill-posed problems in computational vision,” J. Am. Stat. Assoc., vol. 82, no. 397, pp. 76–89, Mar. 1987. [33] M. Kass, A. Withkin, and D. Terzopoulos, “Snakes: Active contour models,” Int. J. Comput. Vis., vol. 1, no. 4, pp. 321–331, Jan. 1988. [34] J. A. Sethian, Level Set Methods. Cambridge, U.K.: Cambridge Univ. Press, 1996. [35] J. F. Aujol, G. Aubert, and L. Blanc-Féraud, “Wavelet-based level set evolution for classification of textured images,” IEEE Trans. Image Process., vol. 12, no. 12, pp. 1634–1641, Dec. 2003. [36] C. Samson, L. Blanc-Féraud, G. Aubert, and J. Zerubia, “A level set method for image classification,” Int. J. Comput. Vis., vol. 40, no. 3, pp. 187–197, 2000. [37] S. Jehan-Besson, M. Barlaud, and G. Aubert, “Image segmentation using active contours: Calculus of variations or shape gradients?” SIAM J. Appl. Math., vol. 63, no. 6, pp. 2128–2154, 2003. [38] H. K. Zhao, T. Chan, B. Merriman, and S. Osher, “A variational level set approach to multiphase motion,” J. Comput. Phys., vol. 127, no. 1, pp. 179–195, Aug. 1996. [39] I. Karoui, R. Fablet, J. M. Boucher, and J. M. Augustin, “Separabilitybased Kullback divergence weighting and filter selection for texture classification and segmentation,” in Proc. 18th ICPR, Hong Kong, Aug. 2006. [40] I. Karoui, R. Fablet, J. M. Boucher, W. Pieczynski, and J. M. Augustin, “Fusion of textural statistics using a similarity measure: Application to texture recognition and segmentation,” Pattern Anal. Appl., vol. 11, no. 3/4, pp. 425–434, Sep. 2008. Computer Science. [41] A. Ehrhold, D. Hamon, and B. Guillaumont, “The Rebent monitoring network, a spatially integrated, acoustic approach to surveying nearshore macrobenthic habitats: Application to the bay of Concarneau (South Brittany, France),” ICES J. Mar. Sci., vol. 63, no. 9, pp. 1604– 1615, Jan. 2006. [42] W. Pieczynski, “Statistical image segmentation,” Mach. Graph. Vis., vol. 1, no. 2, pp. 261–268, 1992. 545 546 547 548 549 550 551 552 553 554 555 556 557 558 559 560 561 562 563 564 565 566 567 568 569 570 571 572 573 574 575 576 577 578 579 580 581 582 583 584 Imen Karoui received the M.Sc. degree in telecommunications from École Polytechnique de Tunis, Tunisia, in 2002, the D.E.A. degree in STIR (signal, télécommunications, images, radar) from the Université de Rennes 1, Rennes, France, in 2003, and the Ph.D. degree in signal processing and telecommunications from the Department of Signal and Communications, École Nationale Supérieure des Télécommunications de Bretagne (Telecom Bretagne), Brest Cedex 3, France, in 2007. She is currently a Postdoctor with Telecom Bretagne. Her research interests include texture analysis and segmentation with application to acoustic imaging and fisheries. 585 586 587 588 589 590 591 592 593 594 595 596 597 Ronan Fablet receive the M.Sc. degree in signal processing and telecommunications from École Nationale Supérieure de l’Aéronautique et de l’Espace, Toulouse Cedex 4, France, in 1997 and the Ph.D. degree in signal processing and telecommunications from the University of Rennes 1, Rennes, France, in 2001. In 2002, he was an INRIA Postdoctoral Fellow with Brown University, Providence, RI. From 2003 to 2007, he had a full-time research position with the Institut Français de Recherche et d’Exploitation de la Mer, Brest, France, in the field of biological oceanography. Since 2008, he has been an Assistant Professor with the Department of Signal and Communications, École Nationale Supérieure des Télécommunications de Bretagne, Brest Cedex 3, France. His main interests include sonar and radar imaging particularly applied to biological and physical oceanography and signal and image processing applied to biological archives. 598 599 600 601 602 603 604 605 606 607 608 609 610 611 612 613 614 IE E Pr E oo f 468 469 470 471 472 473 474 475 476 477 478 479 480 481 482 483 484 485 486 487 488 489 490 491 492 493 494 495 496 497 498 499 500 501 502 503 504 505 506 507 508 509 510 511 512 513 514 515 516 517 518 519 520 521 522 523 524 525 526 527 528 529 530 531 532 533 534 535 536 537 538 539 540 541 542 543 544 IEEE TRANSACTIONS ON GEOSCIENCE AND REMOTE SENSING, VOL. 46, NO. 12, DECEMBER 2008 KAROUI et al.: SEABED SEGMENTATION USING OPTIMIZED STATISTICS OF SONAR TEXTURES Jean-Marc Boucher (M’83) was born in 1952. He received the M.Sc. degree in telecommunications from École Nationale Supérieure des Télécommunications, Paris, France, in 1975 and the Habilitation à Diriger des Recherches degree from the University of Rennes 1, Rennes, France, in 1995. He is currently a Professor with the Department of Signal and Communications, École Nationale Supérieure des Télécommunications de Bretagne, Brest Cedex 3, France, where he is also an Education Deputy Director. He is also the Deputy Manager of a National Scientific Research Center Laboratory (Institut Telecom-Telecom Bretagne-UMR CNRS 3192 lab-STICC, Université Européenne de BretagneFrance). His current researches are related to statistical signal analysis, including estimation theory, Markov models, blind deconvolution, wavelets, and multiscale-image analysis. These methods are applied to radar and sonar images, seismic signals, and electrocardiographic signals. He published about a hundred technical articles in these areas in international journals and conferences. Jean-Marie Augustin received the Master degree in signal processing from the University of Rennes 1, Rennes, France, in 1982. Since 1984, he has been with the Institut Français de Recherche et d’Exploitation de la Mer, Brest, France, the French public institute for marine research, where he is currently a Senior Engineer with the Department of Acoustics and Seismics. His main interests include software development for sonar seafloor mapping and backscatter reflectivity analysis. IE E Pr E oo f 615 616 617 618 619 620 621 622 623 624 625 626 627 628 629 630 631 632 633 11 634 635 636 637 638 639 640 641 642 643 644 IEEE TRANSACTIONS ON GEOSCIENCE AND REMOTE SENSING, VOL. 46, NO. 12, DECEMBER 2008 2 Seabed Segmentation Using Optimized Statistics of Sonar Textures 3 Imen Karoui, Ronan Fablet, Jean-Marc Boucher, Member, IEEE, and Jean-Marie Augustin 1 Abstract—In this paper, we propose and compare two supervised algorithms for the segmentation of textured sonar images, with respect to seafloor types. We characterize seafloors by a set of empirical distributions estimated on texture responses to a set of different filters. Moreover, we introduce a novel similarity measure between sonar textures in this feature space. Our similarity measure is defined as a weighted sum of Kullback–Leibler divergences between texture features. The weight setting is twofold. First, each filter is weighted according to its discrimination power: The computation of these weights are issued from a margin maximization criterion. Second, an additional weight, evaluated as an angular distance between the incidence angles of the compared texture samples, is considered to take into account sonar-image acquisition process that leads to a variability of the backscattered value and of the texture aspect with the incidence-angle range. A Bayesian framework is used in the first algorithm where the conditional likelihoods are expressed using the proposed similarity measure between local pixel statistics and the seafloor prototype statistics. The second method is based on a variational framework as the minimization of a region-based functional that involves the similarity between global-region texture-based statistics and the predefined prototypes. IE E Pr E oo f 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 1 Fig. 1. Typical sidscan sonar image (Rebent, IFREMER). 26 Index Terms—Active regions, angular backscattering, feature 27 selection, level sets, maximum marginal probability (MMP), seg28 mentation, sonar images, texture. I. I NTRODUCTION 29 COUSTIC remote sensing, such as high-resolution multibeam and sidescan sonars, provides new means for in situ observation of the seabed. The characterization of these high33 resolution sonar images is important for a number of practical 34 applications such as marine geology, commercial fishing, off35 shore oil prospecting, and drilling [1]–[4]. 36 The segmentation and the classification of sonar images with 37 respect to seafloor types (rocks, mud, sand, . . .) is the key 38 goal behind the analysis of these acoustic images. This task, 39 however, raises two major difficulties. The first task is to deal 40 with texture in these images. Previous methods are generally 41 based only on backscattered (BS) intensity models, and several 42 parametric families of probabilistic distribution functions have 30 31 32 A Manuscript received January 4, 2008; revised May 20, 2008 and August 11, 2008. I. Karoui, R. Fablet, and J.-M. Boucher are with the Department of Signal and Communications, École Nationale Supérieure des Télécommunications de Bretagne, 29238 Brest Cedex 3, France (e-mail: [email protected]). J.-M. Augustin is with the Department of Acoustics and Seismics, Institut Français de Recherche et d’Exploitation de la Mer, 29280 Brest, France. Color versions of one or more of the figures in this paper are available online at http://ieeexplore.ieee.org. Digital Object Identifier 10.1109/TGRS.2008.2006362 Fig. 2. BS evolution with the incidence angles for the three seafloor types of Fig. 1: sand, mud, and sand ripples. been suggested (Rayleigh distribution, K distribution, Weibull 43 distribution, etc.) [3]–[8]. These first-order statistics are not 44 sufficient when high-resolution sonar images involve textures, 45 which is the case of most sonar images (Fig. 1). 46 The other important issue arising in seabed texture charac- 47 terization is a built-in feature of sonar observation: The value 48 of BS measure depends both on the seafloor type and on the 49 incident angle of the reflected acoustic signal, ranging typically 50 from −85◦ to +85◦ . Fig. 2 shows the BS evolution for three 51 different seafloor types: sand, sand ripples, and mud. In addition 52 to the BS variability within incidence angles, seafloor textures 53 0196-2892/$25.00 © 2008 IEEE 2 IEEE TRANSACTIONS ON GEOSCIENCE AND REMOTE SENSING, VOL. 46, NO. 12, DECEMBER 2008 Fig. 3. Sonar image composed of sand ripples. (Left) [5◦ , 40◦ ]. (Right) [80◦ , 85◦ ]. 55 56 57 for two angular sectors. (Left) [5◦ , 40◦ ]. segmentation schemes. We first state the segmentation issue 93 as a Bayesian pixel-based labeling according to local texture 94 features. The second approach relies on a region-level vari- 95 ational framework, which resorts to a level-set minimization 96 of an energy criterion involving global-region-based seafloor 97 statistics. 98 This paper is organized as follows. The seafloor similarity 99 measure is introduced in Section II. The Bayesian segmentation 100 method is detailed in Section III. The region-based segmen- 101 tation algorithm is described in Section IV. Experiments are 102 reported and discussed in Section V, and conclusions are drawn 103 in Section VI. 104 IE E Pr E oo f are dependent on the incidence angles. Figs. 3 and 4 show a sonar image composed by sand ripples and rocks, respectively, for two angular sectors: [80◦ , 85◦ ] and [5◦ , 40◦ ]. The texture of sand ripples shows a loss of contrast in the specular sec58 tor [5◦ , 40◦ ]: The steep grazing angle reduces the backscatter 59 differences between facing and trailing slopes, while at low 60 incidence angles, much of the variation is lost due to increasing 61 sonar shadow. A similar loss of contrast is observable in the 62 rock samples (Fig. 4). BS behavior according to the incidence 63 angles has been of wide interest for sonar imaging [3], [9]–[12]. 64 Parametric and nonparametric techniques have been proposed 65 to model sonar-image behavior with respect to the incidence 66 angle variations. The effect of the incidence angle on the BS 67 has also been explored as a discriminating feature for seafloor 68 recognition [4], [13], [14]. However, no studies have proposed 69 a method to accurately compensate this phenomena because 70 of the joint dependence of the seafloor types and on the local 71 bathymetry which is generally unknown for sidescan sonar 72 images. To our knowledge, the effects of the incidence angle on 73 textured seabed features have not been addressed for segmen74 tation issues. Only some studies were interested in simulating 75 the behavior of oriented and textured seafloor types [15]–[18]. 76 These methods are mainly based on shape-from-shading [19] 77 and were restricted to synthetic images or to real sonar images 78 involving only one seafloor type. 79 In this paper, we aim at using texture information within 80 sonar seabed images and at developing segmentation algo81 rithms that take into account angular variability of BS and 82 textural features. We propose to characterize seafloors by a 83 wide set of marginal distributions of their filter responses, and 84 we measure seafloor similarities according to a weighted sum 85 of Kullback–Leibler divergences [28] in this feature space. 86 To cope with seafloor angular dependence, we introduce an 87 additional weighting factor, evaluated as an angular distance 88 between the compared texture samples: This angular distance 89 is measured according to Gaussian kernels, whose variance 90 sets the level of the angular variability, depending on textures 91 and seafloor types. The proposed incidence-angle-and-texture92 based similarity measure is exploited to develop two different 54 Fig. 4. Rock texture (Right) [80◦ , 85◦ ]. II. S ONAR -T EXTURE S IMILARITY M EASURE 105 Texture-based segmentation of seabed sonar images gener- 106 ally relies on Haralick [30] parameters or scalar spectral and 107 filter coefficients to model textures [20]–[22]. Recently, in the 108 field of texture analysis, features computed as statistics of local 109 filter responses have been shown to be relevant and discriminant 110 texture descriptors [23]–[27]. Motivated by these studies, we 111 propose to use texture features computed as marginal distrib- 112 utions of a wide set of filter seafloor responses. Each seafloor 113 type denoted by Tk is characterized by a set Qk composed of 114 the marginal distributions of the seafloor with respect to the 115 predefined filters. The following conditions are issued. 116 1) 121 cooccurrence distributions with horizontal and ver- 117 tical displacements denoted by dx and dy , respectively, 118 (dx , dy ) ∈ {0, 1, . . . , 10}. 119 2) 50 distributions of the magnitude of Gabor filter 120 responses, computed for combination of parameters 121 (f0 , σ, θ), where f0 is the radial frequency, σ is 122 the standard√deviation, and θ is the √ orientation, such 123 that f0 ∈ { 2/2k }k=1:6 , σ ∈ {2k 2}k=2:5 , and θ ∈ 124 {0◦ , 25◦ , 45◦ , 90◦ , 135◦ }. 125 3) 48 distributions of the energy of the image wavelet packet 126 coefficient computed for different bands (we used three 127 wavelet types: Haar, Daubechies, and Coiflet). 128 KAROUI et al.: SEABED SEGMENTATION USING OPTIMIZED STATISTICS OF SONAR TEXTURES 3 Gaussian mixture, and the assignment likelihood p(Θj /zs ) is 147 given by 148   −(Θj −θs )2 exp σ2 .  j p(Θj /zs ) =  (1) −(Θi −θs )2 J exp 2 i=1 σ i Hence, the assignment of a given region R to an angular 149 subdomain is given by 150  1 p(Θj /zs ), ds. (2) πR (j) = |R| s IE E Pr E oo f Fig. 5. Cooccurrence distributions computed for parameters (d, θ) = (6, 0◦ ). (Left) [5◦ , 40◦ ]. (Right) [80◦ , 85◦ ]. The angular bounds, defined by the parameters Θj and σj , 151 are set experimentally, and each seafloor type is character- 152 ized by Qk = {Qf,j (Tk )}f =1:F,j=1:J , its filter responses es- 153 timated on the J angular domains. Note that f accounts 154 both for filter types and various associated parameteri- 155 zations: f = 1 → 121 refers to cooccurrence distributions, 156 f = 122 → 171 to Gabor filter, and f = 172 → 219 to 157 wavelet-based distributions. These statistics are computed us- 158 ing Parzen window estimation [29]. Formally, Qf,j (Tk , α) = 159 (1/πTk (j)) Tk p(Θj /zs )gσf (hf (s) − α)ds, where hf is the 160 filter response indexed by f (for a cooccurrence matrix with pa- 161 rameters μ = (d, θ), hf : Ω → [1, N g] × [1, N g], and hf (s) = 162 (I(s), I(s + μ)), where I(s) is the gray value at pixel s and N g 163 is the total gray level number) and gσf is a Gaussian kernel with 164 zero mean and standard deviation σf . 165 We define the similarity measure between a texture sample 166 T with feature set D(T ) = {Df,j (T )}f =1:F,j=1:J and a given 167 seafloor type Tk as follows: 168 J  k  KLΘ w Q , D(T ) = Fig. 6. Cooccurrence distributions computed for parameters (d, θ) = (Left) [5◦ , 40◦ ]. (Right) [80◦ , 85◦ ]. (6◦ , 0). Sonar-texture variability with respect to the incidence angles induces a variability in texture features. Figs. 5 and 6 show the cooccurrence matrices of the images involving sand ripples 132 and rocks shown in Figs. 3 and 4, respectively. A change in 133 cooccurrence distributions can be noticed. In the angular sector 134 [5◦ , 40◦ ], the cooccurrence is bimodal due to the alternation of 135 dark and light values in the image of sand ripples, whereas in 136 the sector [80◦ , 85◦ ], the loss of contrast between dark and light 137 values in the image of sand ripples leads to unimodal distribu138 tion. For cooccurrence distributions related to rock samples, a 139 change in the variance is depicted. 140 To deal with this problem, we propose to define angular 141 subdomains, in which texture characteristics can be regarded 142 as homogeneous. Each angular sector is indexed by j and 143 is characterized by a mean incidence-angle value Θj . For144 mally, this leads in introducing a state variable zs assign145 ing pixel s to a given angular subdomain. (θs , zs ), where 146 θs is the incidence angle of pixel s, is then modeled as a 129 130 131 F j=1 f =1   πT (j)wf2 KL Qkf,j , Df,j (T ) (3) where {wf2 }f =1:F are the feature weights such that 169 F 2 f =1 wf = 1, Df,j (T ) is the distribution related to 170 filter f and estimated on  region T for angular sector j: 171 Df,j (T, α) = (1/πT (j)) T p(Θj /zs )gσf (hf (s) − α)ds, and 172 KL denotes for the Kullback–Leibler divergence [28]; 173 for  two probability distributions Q and D, KL(Q, D) = 174 Q(α)log(Q(α)/D(α))dα. 175 The resulting weighting factors {wf2 }f =1:F are exploited on 176 the one hand for filter selection to keep only the distributions 177 corresponding to the significant weights and, on the other hand, 178 for the definition of an optimized texture-based similarity mea- 179 sure KLΘ w given the selected distributions. In a supervised con- 180 text, the weights are estimated from a training set T composed 181 of N -labeled texture samples: T = {(D(T ), sT )}, where sT is 182 the class of sample T . Formally, {wf2 }f =1:F are issued from the 183 maximization of the global margin defined as follows: 184 MwT = Mw (T ) (4) T ∈T where  d  Θ sT T Mw (T ) = KLΘ w Q , D(T ) − KLw (Q , D(T )) 185 (5) 4 186 187 IEEE TRANSACTIONS ON GEOSCIENCE AND REMOTE SENSING, VOL. 46, NO. 12, DECEMBER 2008 where dT is the nearest class to T different from sT according to the similarity measure KLΘ w  k  (6) dT = arg min KLΘ w Q , D(T ) . k=sT iterative conditional estimation procedure [42]. Using Bayes 227 rule, the posterior distribution is expressed as follows: 228 p(X, Y ) = p(Y /X)p(X) = p(X) The maximization of the margin criterion MwT is carried out using a stochastic gradient method as detailed in our previous 190 work [39]. 188 189 191 III. B AYESIAN S ONAR -I MAGE S EGMENTATION As far as Bayesian image segmentation is concerned, the most popular criteria are the maximum a posteriori (MAP) 194 [31] and the maximum marginal probability (MMP) [32]. It 195 has been shown that the MMP estimation criterion is more 196 appropriate for image segmentation than the MAP criterion 197 [32]. The MAP estimate assigns the same cost to every incorrect 198 segmentation regardless the number of pixels at which the 199 estimated segmentation differs from the true one, whereas the 200 MMP algorithm minimizes the expected value of the number of 201 misclassified pixels. As shown in [32], the MMP procedure is 202 equivalent to maximizing the marginal of the class labels. Let 203 us introduce the following notations: 204 1) S, the image lattice composed of N pixels; 205 2) X = {xs }s∈S , the label random field; 206 3) Y = {ys }s∈S , the random field of the observations, the 207 textural feature in our case. 208 Formally, the MMP scheme is equivalent to the maximi209 zation of x̂MMP = arg max p(xs = k/y). s k∈{1:K} (7) In general, pixel conditional likelihood is computed according to local texture features like Gabor or wavelet coefficients. Here, we aim at using the proposed similarity measure KLΘ w. We associate to each pixel s a set of features D(Ws ) = 214 {Df (Ws )}f =1:F estimated according to a Parzen estimation 215 method [29], within a square window Ws centered at s. The 216 window size that we denote by TW is set by the user according 217 to texture coarseness. In our case, the observation denoted 218 by a random field Y is specified by {ys = D(Ws )}s∈S . The 219 conditional likelihood at each pixel s with respect to class k is 220 then defined from the similarity measure KLΘ w by 210 211 212 213 k −KLΘ w (Q ,D(Ws )) 221 222 223 exp . p(ys /xs = k) = K i −KLΘ w (Q ,D(Ws )) i=1 exp s∈S t∈cs −αc (1−δ(xs ,xt ))+log(p(ys /xs )) . (10) The maximization of local probabilities p(xs = k/y) is car- 229 ried out as follows [32]: We use the Gibbs sampler to gen- 230 erate a discrete-time Markov chain X(t) which converges in 231 distribution to a random field with probability mass functions 232 P (X/Y ). The marginal conditional distributions p(xs = k/y) 233 are then approximated as the fraction of time the Markov chain 234 spends in state k for each class k and each pixel s. The MMP 235 segmentation steps are the following. 236 1) Simulation of Tmax realizations of x1 , x2 , . . . , xTmax of 237 X using Gibbs sampler. 238 2) Using the realizations x1 , x2 , . . . , xTmax , p(xs /y) is esti- 239 mated using the frequency of each realizations 240     δ x1s − k + · · · + δ x1s − k p(xs = k/y) = . Tmax s∈S t∈cs αc (1 − δ(xs , xt )) 3) Choose as xs the class that maximizes p(xs = k/y). 241 IV. V ARIATIONAL S ONAR -I MAGE S EGMENTATION 242 Unlike the Bayesian scheme, the second approach is stated 243 at a region level as the minimization of a constrained energy 244 criterion E({Ωk }k=1:K ) = E1 + E2 , where Ωk is the domain 245 composed of all pixels attributed to the class k, E1 is a texture- 246 based data-driven term, and E2 is a regularization term detailed 247 as follows. 248 A. Functional Terms E ({Ωk }k=1:K ) is evaluated as the log-likelihood of a given 250 partition with respect to texture models. It is evaluated as the 251 sum of the similarities according to the measure KLΘ w between 252 each region Ωk and its corresponding class Tk 253 K E1 ({Ωk }k=1:K ) = k=1  k  KLΘ w Q , D(Ωk ) (11) where Df,j (Ωk ) is the marginal distribution of the image 254 response to the filter indexed by f . For angular domain j, 255 Df,j (Ωk ) is estimated according to Parzen method [29]. 256 E2 is a regularization term, it penalizes the lengths of region 257 contours and is expressed by 258 (9) where δ is the delta function and αc ∈ {αH , αV , αD } are real parameters assigned, respectively, to horizontal, vertical, and 226 diagonal cliques. These parameters are estimated using the 249 1 (8) As a prior PX , we consider a Markov random field associated to an eight-neighborhood system with potential functions given by U2 (x) = 224 225 ≈ exp   IE E Pr E oo f 192 193 p(ys /xs ) s K E2 = k=1 γk |Γk |, γk ∈ ℜ+ (12) where |Γk | is the length of the contour Γk associated to the 259 region Ωk . 260 KAROUI et al.: SEABED SEGMENTATION USING OPTIMIZED STATISTICS OF SONAR TEXTURES 261 B. Computation of the Evolution Equation We solve for the minimization of the functional E using a gradient-descent technique. It relies on the computation of 264 the first derivative of E according to regions {Ωk }k=1:K . The 265 evolution equation of region contours {Γk }k=1:K is then given 266 by the following dynamic scheme [33]: 262 263 ∂Γk (x,t) = Fk (x, t)Nk ∂t Γk (x, 0) = Γ0k (13) 5 where ∂ϕik /∂t, i = 1, 2, 3 are the evolution equation terms associated, respectively, to functionals Ei , i = 1, 2, 3. The derivatives of the energy terms E2 and E3 are directly estimated from level-set functions [35], [36] 287 288 289 290 ∇ϕk ∂ϕ2k (s, t) = γk δα (ϕk )div ∀k ∈ {1 : K} (19) ∂t |∇ϕk | K  ∂ϕ3k (s, t) =−δα (ϕk )λ (Hα (ϕk )−1) ∀k ∈ {1 : K}. ∂t k=1 where Nk is the unit inward normal to Γk at pixel x and at time t and Fk is the velocity field (in our case, Fk = ∇Ek , the 269 derivative of E with respect to Γk ). 270 The explicit implementation of the curve evolution according 271 to the latter dynamic scheme using a difference-approximation 272 scheme cannot deal with topological changes of the moving 273 front. This could be avoided by introducing the level-set method 274 proposed by Osher and Sethian [34]. The basic idea of the 275 method is the implicit representation of the moving interface Γ 276 by a higher dimensional hypersurface ϕ (the level-set function) 277 such that the zero level set of ϕ is actually the set of Γ and 267 268 K γk lim k=1 281 ∀s ∈ Ωk . (14) E2 can be expressed using level-set functions ϕk [36] E2 = α→0  δα (ϕk )|∇ϕk |ds (15) Ω where δα is a regularized delta function    1 1 + cos πs , 2α α δα (s) = 0, if |s| ≤ α if |s| < α. (16) In order to cope with multiclass segmentation and to fulfill the image-partition constraint, we use an additional term E3 given 284 by the following functional [38]: 282 283 λ E3 = 2   Ω 285 286 K k=1 F f =1    k wf2 p(Θj /zs ) KL Qkf,j , Df,j  − Qkf,j /Df,j ∗ gσf (hf (s)) − 1 |∇ϕk | The evolution of the contours {Γk }k=1:K (13) is then equivalent to the evolution of the level-set functions {ϕk }k=1:K [34] ∂ϕk (s, t) = Fk (s, t) |∇ϕk (s, t)| ∂t 280 J ∂ϕ1k (s, t) = − ∂t j=1 IE E Pr E oo f 278 The evolution equation related to E1 is more complex, since it 291 involves computations over the spatial support of each region. 292 To differentiate E1 , we use shape-derivative tools, particularly 293 the Gâteaux derivative theorem given in [37]. As detailed in the 294 Appendix, it leads to 295  Ωinside = {s ∈ Ω/ϕ(s) > 0} Ωoutside = {s ∈ Ω/ϕ(s) < 0} . 279 (20) 2 Hα (ϕk ) − 1 ds, λ ∈ ℜ+ where Hα is a regularized Heaviside function ⎧    ⎨ 12 1 + αs + π1 sin πs , if |s| ≤ α α Hα (s) = 1, if s > α ⎩ 0, if s < −α. As E = E1 + E2 + E3 , we have ∂ϕk ∂ϕ1k ∂ϕ2k ∂ϕ3k (s, t) = (s, t) + (s, t) + (s, t) ∂t ∂t ∂t ∂t (17) (21) where ∗ is the convolution symbol. 296 The evolution equation related to the energy E1 (21) is 297 composed of two terms. 298 J  F 2 k 1) A global term − j=1 f =1 wf p(Θj /zs )KL(Qf,j , 299 Df,j ): This term is always negative or null. It is a contrac- 300 tion force that  reducesthe size of heterogeneous regions. 301 F J 2 k 2) A local term f =1 wf p(Θj /zs )((Qf,j /Df,j ) − 302 j=1 1) ∗ gσf (hf (s)): This term locally compares the features 303 values at each pixel. This term can be positive or neg- 304 ative and aims at readjusting the statistics inside the 305 regions Df,j (Ωk ) to fit to prototype models Qkf,j . The 306 contribution of each descriptor f is weighted by wf2 and 307 p(Θj /zs ), the relative contribution of descriptor f and of 308 angular sector j. 309 The overall evolution equations of the contours {Γk }k=1:K 310 are the following: 311 ∂ϕk (s, t) ∂t = −δα (ϕk ) ⎡ J ×⎣ F wf2 p(Θj /zs ) j=1 f =1 ×  KL(Qkf,j , Df,j )− (18) + γk div ∇ϕk −λ |∇ϕk |  Qkf,j ∗ gσf (hf (s))+1 Df,j  K (Hα (ϕk )−1) k=1 ∀k. (22) We apply these coupled evolution equations until convergence. 312 6 IEEE TRANSACTIONS ON GEOSCIENCE AND REMOTE SENSING, VOL. 46, NO. 12, DECEMBER 2008 Fig. 8. Feature weights. 1 → 121 cooccurrence distributions, 122 → 171 Gabor filter-based features and 172 → 219 wavelet-based distributions. (a) I1. (b) I5. IE E Pr E oo f TABLE I SEGMENTATION ERROR RATES Fig. 7. Test images and their manual segmentation (in black line). 313 V. E XPERIMENTAL R ESULTS AND D ISCUSSION In previous work, we have tested the method on various optic textures (Brodatz textures). The method was compared 316 to other texture-classification methods, and some results are 317 reported in [39] and [40]. Here, we evaluate the proposed 318 seabed-segmentation technique for different real sonar images 319 acquired by a sidescan sonar, as part of a natural seabed320 mapping project (IFREMER, Rebent Project) [41]. A reference 321 interpretation by an expert is available [41]. Fig. 7 shows the 322 set of images on which we carried out the experiments. We 323 superimposed on these images the manual expert segmentation. 324 Image I1 is composed of three seafloor types: rock, mud, and 325 marl ripples [41]; I2 of rock, marl ripples, and mud seafloors; 326 I3 of mud, sand, and marl ripples; I4 of marl and marl ripples; 327 and I5 of sand, sand ripples, and rock. For I2 and I3, the angular 328 variability of the sefloors, particularly marl ripples (for I2) and 329 marl (for I3), is visually clear. 330 For all these images, we first determine the most discriminant 331 features among the initial set of 219 features: We apply the 332 algorithm described in Section II and detailed in [39], and 333 we keep only the feature set such that the cumulative sum 334 of weights exceeds 0.9. Only a small number of features are 335 retained. For example, Fig. 8 shows the plot of feature weights 336 computed for image I1 and I5. For I1, the two cooccurrence 337 matrices account for more than 90% of the total weight sum; 338 these cooccurrence are computed for parameters (dx , dy ) ∈ 339 {(1, 4), (2, 1)}. For I5, the cooccurrence distributions com340 puted for parameters (dx , dy ) ∈ {(2, 1), (2, 2)} are selected. 341 For sonar images, we noticed that cooccurrence matrices are the 314 315 most selected features. In previous work on Brodatz textures 342 [39], [40], we remarked that Gabor and wavelet filters were 343 selected for oriented textures, whereas cooccurrence distribu- 344 tions, which, in addition to the detection of texture structures, 345 detect the intensity values change, are selected in the case of 346 texture having different intensity values and texture with regular 347 motifs. 348 For the five test images, three segmentation algorithms are 349 compared. 350 1) The maximum-likelihood (ML) segmentation denoted by 351 ML. This method consists in maximizing at each pixel the 352 conditional probability p(ys /xs ) given by (8) 353 xˆs = arg max p(ys /xs = k). k (23) Several analysis-window sizes are compared: TW ∈ {7 × 354 7, 17 × 17, 33 × 33}. 355 2) The MMP segmentation described in Section III, applied 356 for several analysis-window sizes: TW ∈ {7 × 7, 17 × 357 17, 33 × 33}. 358 3) The region-based variational segmentation described in 359 Section IV that we denote by V ar. 360 Table I summarizes the different error classification rates for all 361 segmentations. 362 All segmentation methods give quite good results according 363 to the mean classification error rates. MMP and variational 364 approaches are more efficient than the ML-based segmentation 365 because they take into account the spatial dependence between 366 pixels. The difference between these later approaches (MMP 367 and variational methods) mainly lies in the accuracy of the 368 localization of region boundaries and in the dependence of 369 MMP-based segmentation on the window size TW . In fact, 370 for the MMP segmentation, the use of small neighborhood 371 (TW = 7 × 7) leads to more accurate region frontiers but 372 small misclassified patches appear because of the neighborhood 373 KAROUI et al.: SEABED SEGMENTATION USING OPTIMIZED STATISTICS OF SONAR TEXTURES MMP segmentation of I1 using TW = 7 × 7, τ = 12.5%. Fig. 10. MMP segmentation of I1 using TW = 17 × 17, τ = 9.2%. inefficiency for texture characterization. Conversely, the use of large window sizes (TW = 33 × 33) resorts to a lack of 376 accuracy in the localization of the boundaries of the seabed 377 regions because texture features extracted for pixels close to 378 region boundaries involve a mixture of texture characteris379 tics. The variational region-based approach does not need the 380 choice of an analysis window and operates globally on region 381 composed of pixels belonging to the same class. It resorts 382 to a tradeoff between segmentation accuracy and region ho383 mogeneity. Figs. 9–13 show examples of the dependence of 384 MMP segmentation on the sizes of analysis window and on the 385 robustness of the variational approach. In Fig. 14, we plot the 386 mean error rate for different window sizes. 387 We note that MMP and variational segmentations give sim388 ilar results when the analysis-window size is well chosen for 389 instance, TW = 17 × 17 for image I2 (Figs. 15 and 16), but 390 its performance depends a lot on the choice of this parameter. 391 This method can however be appropriate when the aim of the 392 segmentation is to detect texture regions and without seeking 393 for accurate boundaries. 394 The variational approach is also interesting because it is 395 much faster than the MMP segmentation. Being deterministic, 396 the variational approach can be very fast particularly if we use 397 appropriate initialization such as an initial segmentation based 398 on the ML criterion. Whereas the MMP segmentation needs 374 375 Fig. 11. MMP segmentation of I1 using TW = 33 × 33, τ = 13%. IE E Pr E oo f Fig. 9. 7 Fig. 12. Region-based segmentation of I1, τ = 7.5%. a large number of iterations, each iteration is also complex 399 and involves Gibbs sampling. For our implementations, the 400 convergence time of a variational image typically corresponds 401 to one iteration of the MMP algorithm. 402 To stress the interest of taking into account texture variability 403 with respect to incidence angles, additional segmentation re- 404 sults are reported for image I2 and I3 for which the seafloor 405 texture variability is clearer. In Table II, we summarize the 406 segmentation error rates for the segmentation with (J = 3) and 407 with no angular weighting, i.e., using only one angular domain 408 J = 1, for ML, MMP, and variational methods. In Figs. 17 and 409 18 are reported the associated segmentation results. It can be 410 noticed that a classical segmentation (without taking into ac- 411 count texture variability within incidence angles) cannot distin- 412 guish between visually similar seafloors (mud and marl ripples 413 for I2 and marl and sand for I3) near the specular domain. 414 VI. C ONCLUSION 415 We proposed two segmentation algorithms for sonar-image 416 segmentation: a Bayesian algorithm using local statistics and a 417 region-based variational algorithm, both based on a novel sim- 418 ilarity measure between seafloor-type images according to the 419 statistics of their responses to a large set of filters. This similar- 420 ity measure is expressed as a weighted sum of Kullback–Leibler 421 8 IEEE TRANSACTIONS ON GEOSCIENCE AND REMOTE SENSING, VOL. 46, NO. 12, DECEMBER 2008 IE E Pr E oo f Fig. 13. Segmentations of I4. (a) MMP: TW = 7, τ = 4%. (b) MMP: TW = 33, τ = 6%. (c) Variational method: τ = 3%. Fig. 16. Variational-based segmentation of I2, τ = 3%. Fig. 14. Mean segmentation error rate for several window sizes. TABLE II SEGMENTATION ERROR RATES Fig. 15. MMP-based segmentation of I2 TW = 17, τ = 4%. divergences between individual seafloor filter-response statistics. The resulting weighting factors are exploited on the one hand for filter selection and, on the other hand, for taking into account the incidence angular dependence of seafloor tex426 tures. The conclusion is the cooccurrence matrices outperform 427 the other features for our sonar images. The results show 428 that the performance of the Bayesian approach depends on 429 the size of the analysis window. For pixel-based segmenta430 tion (ML and MMP segmentations), the size must be tuned 431 according to the coarseness of given textures. The results 432 also stress the suitability of the region-based approach as 433 compared to the Bayesian pixel-based scheme for texture seg434 mentation and demonstrate the interest of taking into account 422 423 424 425 Fig. 17. I2 region-based segmentation with and without angular weighting. (a) J = 1, τ = 23.5%. (b) J = 3, τ = 3%. Fig. 18. I3 MMP-based segmentation with and without angular weighting. (a) J = 1, τ = 16%. (b) J = 3, τ = 8%. the angular backscatter variabilities to discriminate between 435 seafloor types particularly near the specular sector. 436 KAROUI et al.: SEABED SEGMENTATION USING OPTIMIZED STATISTICS OF SONAR TEXTURES A PPENDIX E VOLUTION E QUATION C OMPUTATION 437 438 441 Using the shape-derivative tools, we want to differentiate the functional  k  F (Ωk ) = KLΘ w Q , D(Ωk ) which can be written as follows: J F πΩk (j)wf2 j=1 f =1 442 443 444  Qkf,j (α) log Rf Qkf,j (α) Df,j (Ωk , α)  Then, we have J F wf2 πΩk (j) j=1 f =1 J F Finally, we get 457 ) dHjk (Ωk , V   = p(Θj /zs )  Γk Qkf,j (α) Qf (α) gσf (hf (s) − α) − G1 (f, α, Ωk ) πΩk (j) )  Hjk . Γk + πΩk (j) Hjk dα  )dα. dHjk (Ωk , V  · Nk )da(s) × (V   )dα dHjk (Ωk , V Rf 1 = πΩk (j) Rf    · Nk )da(s) × (V    Qkf,j (α) 1 = gσ (hf (s)−α) Qf (α) p(Θj /zs ) πΩk (j) Df,j (Ωk , α) f Rf  wf2 dπΩk (j)(Ωk , V j=1 f =1 Qkf,j ∂Hjk = . ∂πΩk (j) πΩk (j) Γk The Gâteaux derivative of F (Ωk ) in the direction of a vector  is given by field V )= dF (Ωk , V Qkf,j ∂Hjk =− ∂G1 G1 According to the previous theorem, we get 456   · Nk )da(s).  ) = − p(Θj /zs )gσ (hf (s) − α) (V dG1 (Ωk , V f p(Θj /zs )gσf (hf (s) − α) ds. F (Ωk ) = 446 dα. 2) Ωk 447  Let us introduce the following notations. 1)   Qkf,j (α) k k Hj = Qf,j (α) log . Df,j (Ωk , α) G1 = 445  455 Using the chain rule, we get   ∂Hjk ∂Hjk  )=  )+ ) dHjk (Ωk , V dG1(Ωk , V dπΩk(j)(Ωk , V ∂G1 ∂πΩk (j) IE E Pr E oo f 439 440 9  Γk   Qkf,j ∗ gσf (hf (s)) − 1 p(Θj /zs ) Df,j (Ωk )  · Nk )da(s). × (V (24) Finally, according to (24) 458 Rf 448 Theorem: See [37]. 449  The Gâteaux derivative 450 Ω k(s, Ω)ds is given by )= dK(Ω, V  ) dF (Ωk , V of a functional of the type K(Ω) =  )ds − ksh (s, Ω, V  N  )da(s) k(s, Ω)(V Γ Ω J =− F p(Θj /zs )wf2 j=1 f =1       k  Qkf,j × KL Qf,j , Df,j − −1 ∗ gσf (hf (s)) Df,j (Ωk ) Γk  ) is where Ω denoted a region and Γ its boundary; ksh (s, Ω, V the shape derivative of k(s, Ω). 453 In our case   ) = − p(Θj /zs )(V  · Nk )da(s) dπΩk (j)(Ωk , V 451 452 Γk 454 Hjk can be written as Hjk = Qkf,j (α) log  πΩk (j)Qkf,j (α) G1 (Ωk , f, j, α)  .    · Nk da(s). × V R EFERENCES 459 [1] P. Brehmer, F. Gerlotto, J. Guillard, F. Sanguinède, Y. Guénnegan, and D. Buestel, “New applications of hydroacoustic methods for monitoring shallow water aquatic ecosystems: The case of mussel culture grounds,” Aquat. Living Resour., vol. 16, no. 3, pp. 333–338, Jul. 2003. [2] M. Mignotte, C. Collet, P. Perez, and P. Bouthemy, “Hybrid genetic optimization and statistical model based approach for the classification of shadow shapes in sonar imagery,” IEEE Trans. Pattern Anal. Mach. Intell., vol. 22, no. 2, pp. 129–141, Feb. 2000. 460 461 462 463 464 465 466 467 10 [3] L. Hellequin, J. M. Boucher, and X. Lurton, “Processing of highfrequency multibeam echo sounder data for seafloor characterization,” IEEE J. Ocean. Eng., vol. 28, no. 1, pp. 78–89, Jan. 2003. [4] G. Le Chenadec and J. M. Boucher, “Sonar image segmentation using the angular dependence of backscattering distributions,” in Proc. Oceans—Europe, Brest, France, 2005, vol. 1, pp. 147–152. [5] G. Le Chenadec, J.-M. Boucher, and X. Lurton, “Angular dependence of K-distributed sonar data,” IEEE Trans. Geosci. Remote Sens., vol. 45, no. 5, pp. 1224–1236, May 2007. [6] E. Jakeman, “Non-Gaussian models for the statistics of scattered waves,” Adv. Phys., vol. 37, no. 5, pp. 471–529, Oct. 1988. [7] A. P. Lyons and D. A. Abraham, “Statistical characterization of highfrequency shallow-water seafloor backscatter,” J. Acoust. Soc. Amer., vol. 106, no. 3, pp. 1307–1315, Sep. 1999. [8] D. A. Abraham and A. P. Lyons, “Novel physical interpretations of k-distributed reverberation,” IEEE J. Ocean. Eng., vol. 27, no. 4, pp. 800– 813, Oct. 2002. [9] N. C. Mitchell, “Processing and analysis of Simrad multibeam sonar data,” Mar. Geophys. Res., vol. 18, no. 6, pp. 729–739, Dec. 1996. [10] G. Le Chenadec, J. M. Boucher, X. Lurton, and J. M. Augustin, “Angular dependence of statistical distributions for backscattered signals: Modeling and application to multibeam echosounder data,” in Proc. Oceans, Sep. 2003, vol. 2, pp. 897–903. [11] D. R. Jackson, D. P. Winbrenner, and A. Ishimaru, “Application of the composite roughness model to high-frequency bottom backscattering,” J. Acoust. Soc. Amer., vol. 79, no. 5, pp. 1410–1422, May 1986. [12] C. de Moustier and D. Alexandrou, “Angular dependence of 12-kHz seafloor acoustic backscatter,” J. Acoust. Soc. Amer., vol. 90, no. 1, pp. 522–531, Jul. 1991. [13] J. E. Hughes-Clarke, “Toward remote seafloor classification using the angular response of acoustic backscatter: A case study from multiple overlapping GLORIA data,” IEEE J. Ocean. Eng., vol. 19, no. 1, pp. 112– 127, Jan. 1994. [14] E. Pouliquen and X. Lurton, “Automated sea-bed classification system for echo-sounders,” in Proc. Oceans, Oct. 1992, vol. 1, pp. 317–321. [15] D. Langer and M. Herbert, “Building qualitative elevation maps from side scan sonar data for autonomous underwater navigation,” in Proc. IEEE Int. Conf. Robot. Autom., Apr. 1991, vol. 3, pp. 2478–2483. [16] R. Li and S. Pai, “Improvement of bathymetric data bases by shape from shading technique using side-scan sonar images,” in Proc. OCEANS, Oct. 1991, vol. 1, pp. 320–324. [17] E. Dura, J. Bell, and D. Lane, “Reconstruction of textured seafloors from side-scan sonar images,” Proc. Inst. Elect. Eng.—Radar Sonar Navig., vol. 151, no. 2, pp. 114–125, Apr. 2004. [18] A. E. Johnson and M. Hebert, “Seafloor map generation for autonomous underwater vehicle navigation,” Auton. Robots, vol. 3, no. 2, pp. 145–168, Jun. 1996. [19] R. Zhang, P. Tsai, J. E. Cryer, and M. Shah, “Shape from shading: A survey,” IEEE Trans. Pattern Anal. Mach. Intell., vol. 21, no. 8, pp. 690– 705, Aug. 1999. [20] T. H. Eggen, “Acoustic sediment classification experiment by means of multibeam echosounders,” Acoust. Classification Mapping Seabed, vol. 15, no. 2, pp. 437–444, 1993. [21] L. M. Linnett, S. J. Clarke, and D. R. Carmichael, “The analysis of sidescan sonar images for seabed types and objects,” in Proc. 2nd Eur. Conf. Underwater Acoust., 1994, vol. 22, pp. 187–194. [22] M. Lianantonakis and Y. R. Petillot, “Sidescan sonar segmentation using active contours and level set method,” in Proc. OCEANS—Europe, Brest, France, 2005, vol. 1, pp. 719–724. [23] L. Xiuwen and W. DeLiang, “Texture classification using spectral histograms,” IEEE Trans. Image Process., vol. 12, no. 6, pp. 661–670, Jun. 2003. [24] O. G. Cula and K. Dana, “3D texture recognition using bidirectional feature histograms,” Int. J. Comput. Vis., vol. 59, no. 1, pp. 33–60, Aug. 2004. [25] R. Fablet and P. Bouthemy, “Motion recognition using nonparametric image motion models estimated from temporal and multiscale co-occurrence statistics,” IEEE Trans. Pattern Anal. Mach. Intell., vol. 25, no. 12, pp. 1619–1624, Dec. 2003. [26] P. Nammalwar, O. Ghita, and P. F. Whelan, “Integration of feature distributions for colour texture segmentation,” in Proc. Int. Conf. Pattern Recog., Aug. 2005, vol. 1, pp. 716–719. [27] Q. Xu, J. Yang, and S. Ding, “Texture Segmentation using LBP embedded region competition,” Electron. Lett. Comput. Vis. Image Anal., vol. 5, no. 1, pp. 41–47, 2004. [28] S. Kullback, “On information and sufficiency,” in The Annals of Mathematical Statistics, vol. 22. New York: Wiley, 1951. [29] E. Parzen, “On estimation of a probability density function and mode,” Ann. Math. Stat., vol. 33, no. 3, pp. 1065–1076, Sep. 1962. [30] R. M. Haralick, “Statistical and structural approaches to texture,” Proc. IEEE, vol. 67, no. 5, pp. 786–804, May 1979. [31] S. Geman and G. Geman, “Stochastic relaxation, Gibbs distributions, and the Bayesian restoration of images,” IEEE Trans. Pattern Anal. Mach. Intell., vol. PAMI-6, no. 6, pp. 721–741, Nov. 1984. [32] J. Marroquin, S. Mitter, and T. Poggio, “Probabilistic solution of ill-posed problems in computational vision,” J. Am. Stat. Assoc., vol. 82, no. 397, pp. 76–89, Mar. 1987. [33] M. Kass, A. Withkin, and D. Terzopoulos, “Snakes: Active contour models,” Int. J. Comput. Vis., vol. 1, no. 4, pp. 321–331, Jan. 1988. [34] J. A. Sethian, Level Set Methods. Cambridge, U.K.: Cambridge Univ. Press, 1996. [35] J. F. Aujol, G. Aubert, and L. Blanc-Féraud, “Wavelet-based level set evolution for classification of textured images,” IEEE Trans. Image Process., vol. 12, no. 12, pp. 1634–1641, Dec. 2003. [36] C. Samson, L. Blanc-Féraud, G. Aubert, and J. Zerubia, “A level set method for image classification,” Int. J. Comput. Vis., vol. 40, no. 3, pp. 187–197, 2000. [37] S. Jehan-Besson, M. Barlaud, and G. Aubert, “Image segmentation using active contours: Calculus of variations or shape gradients?” SIAM J. Appl. Math., vol. 63, no. 6, pp. 2128–2154, 2003. [38] H. K. Zhao, T. Chan, B. Merriman, and S. Osher, “A variational level set approach to multiphase motion,” J. Comput. Phys., vol. 127, no. 1, pp. 179–195, Aug. 1996. [39] I. Karoui, R. Fablet, J. M. Boucher, and J. M. Augustin, “Separabilitybased Kullback divergence weighting and filter selection for texture classification and segmentation,” in Proc. 18th ICPR, Hong Kong, Aug. 2006. [40] I. Karoui, R. Fablet, J. M. Boucher, W. Pieczynski, and J. M. Augustin, “Fusion of textural statistics using a similarity measure: Application to texture recognition and segmentation,” Pattern Anal. Appl., vol. 11, no. 3/4, pp. 425–434, Sep. 2008. Computer Science. [41] A. Ehrhold, D. Hamon, and B. Guillaumont, “The Rebent monitoring network, a spatially integrated, acoustic approach to surveying nearshore macrobenthic habitats: Application to the bay of Concarneau (South Brittany, France),” ICES J. Mar. Sci., vol. 63, no. 9, pp. 1604– 1615, Jan. 2006. [42] W. Pieczynski, “Statistical image segmentation,” Mach. Graph. Vis., vol. 1, no. 2, pp. 261–268, 1992. 545 546 547 548 549 550 551 552 553 554 555 556 557 558 559 560 561 562 563 564 565 566 567 568 569 570 571 572 573 574 575 576 577 578 579 580 581 582 583 584 Imen Karoui received the M.Sc. degree in telecommunications from École Polytechnique de Tunis, Tunisia, in 2002, the D.E.A. degree in STIR (signal, télécommunications, images, radar) from the Université de Rennes 1, Rennes, France, in 2003, and the Ph.D. degree in signal processing and telecommunications from the Department of Signal and Communications, École Nationale Supérieure des Télécommunications de Bretagne (Telecom Bretagne), Brest Cedex 3, France, in 2007. She is currently a Postdoctor with Telecom Bretagne. Her research interests include texture analysis and segmentation with application to acoustic imaging and fisheries. 585 586 587 588 589 590 591 592 593 594 595 596 597 Ronan Fablet receive the M.Sc. degree in signal processing and telecommunications from École Nationale Supérieure de l’Aéronautique et de l’Espace, Toulouse Cedex 4, France, in 1997 and the Ph.D. degree in signal processing and telecommunications from the University of Rennes 1, Rennes, France, in 2001. In 2002, he was an INRIA Postdoctoral Fellow with Brown University, Providence, RI. From 2003 to 2007, he had a full-time research position with the Institut Français de Recherche et d’Exploitation de la Mer, Brest, France, in the field of biological oceanography. Since 2008, he has been an Assistant Professor with the Department of Signal and Communications, École Nationale Supérieure des Télécommunications de Bretagne, Brest Cedex 3, France. His main interests include sonar and radar imaging particularly applied to biological and physical oceanography and signal and image processing applied to biological archives. 598 599 600 601 602 603 604 605 606 607 608 609 610 611 612 613 614 IE E Pr E oo f 468 469 470 471 472 473 474 475 476 477 478 479 480 481 482 483 484 485 486 487 488 489 490 491 492 493 494 495 496 497 498 499 500 501 502 503 504 505 506 507 508 509 510 511 512 513 514 515 516 517 518 519 520 521 522 523 524 525 526 527 528 529 530 531 532 533 534 535 536 537 538 539 540 541 542 543 544 IEEE TRANSACTIONS ON GEOSCIENCE AND REMOTE SENSING, VOL. 46, NO. 12, DECEMBER 2008 KAROUI et al.: SEABED SEGMENTATION USING OPTIMIZED STATISTICS OF SONAR TEXTURES Jean-Marc Boucher (M’83) was born in 1952. He received the M.Sc. degree in telecommunications from École Nationale Supérieure des Télécommunications, Paris, France, in 1975 and the Habilitation à Diriger des Recherches degree from the University of Rennes 1, Rennes, France, in 1995. He is currently a Professor with the Department of Signal and Communications, École Nationale Supérieure des Télécommunications de Bretagne, Brest Cedex 3, France, where he is also an Education Deputy Director. He is also the Deputy Manager of a National Scientific Research Center Laboratory (Institut Telecom-Telecom Bretagne-UMR CNRS 3192 lab-STICC, Université Européenne de BretagneFrance). His current researches are related to statistical signal analysis, including estimation theory, Markov models, blind deconvolution, wavelets, and multiscale-image analysis. These methods are applied to radar and sonar images, seismic signals, and electrocardiographic signals. He published about a hundred technical articles in these areas in international journals and conferences. Jean-Marie Augustin received the Master degree in signal processing from the University of Rennes 1, Rennes, France, in 1982. Since 1984, he has been with the Institut Français de Recherche et d’Exploitation de la Mer, Brest, France, the French public institute for marine research, where he is currently a Senior Engineer with the Department of Acoustics and Seismics. His main interests include software development for sonar seafloor mapping and backscatter reflectivity analysis. IE E Pr E oo f 615 616 617 618 619 620 621 622 623 624 625 626 627 628 629 630 631 632 633 11 634 635 636 637 638 639 640 641 642 643 644