Clustering stability for automated color image segmentation

Accepted Manuscript

Clustering stability for automated color image segmentation Ariel E. Baya, G. Larese, Rafael Nam´ıas ´ Monica ´ PII: DOI: Reference:

S0957-4174(17)30393-7 10.1016/j.eswa.2017.05.064 ESWA 11356

To appear in:

Expert Systems With Applications

Received date: Revised date: Accepted date:

28 October 2016 2 May 2017 27 May 2017

Please cite this article as: Ariel E. Baya, G. Larese, Rafael Nam´ıas, Clustering stabil´ Monica ´ ity for automated color image segmentation, Expert Systems With Applications (2017), doi: 10.1016/j.eswa.2017.05.064

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

ACCEPTED MANUSCRIPT

Highlights • We propose a clustering validation method (CS) to automatically segment

CR IP T

images. • We adapt CS to detect the best settings for color-texture feature extraction.

• Our method provides a score profile to compare alternative segmentation solutions.

AN US

• We test our procedure on textures, color natural images and 3D MRI data.

• Our method outperforms results obtained with other clustering validation

AC

CE

PT

ED

M

indexes

1

ACCEPTED MANUSCRIPT

Clustering stability for automated color image segmentation

a CIFASIS,

2

3

CR IP T

Ariel E. Bay´ aa,∗, M´onica G. Laresea , Rafael Nam´ıasa

1

French Argentine International Center for Information and Systems Sciences, CONICET-UNR (Argentina) Bv. 27 de Febrero 210 Bis, 2000 Rosario, Argentina

Abstract

4 5 6

7

8

validation is rarely used for this purpose. In this work we adapt a clustering

9

validation method, Clustering Stability (CS), to automatically segment images.

10

CS is not limited by image dimensionality nor by the clustering algorithm. We

11

show clustering and validation acting together as a data-driven process able to

12

find the optimum number of partitions according to our proposed color-texture

13

feature representation. We also describe how to adapt CS to detect the best

14

settings required for feature extraction. The segmentation solutions found by

15

our method are supported by a stability score named STI, which provides an

16

objective quantifiable metric to obtain the final segmentation results. Further-

17

more, the STI allows to compare multiple alternative solutions and select the

18

most appropriate according to the index meaning. We successfully test our

19

procedure on texture and natural images, and 3D MRI data.

20

Keywords: Clustering validation, Image segmentation, Clustering stability

21

CE

PT

ED

M

AN US

Clustering is a well-established technique for segmentation. However, clustering

1. Introduction

22

23

beled data needs to be studied. Among the techniques comprising it, clustering

24

algorithms are one of the most used methods. There are countless clustering

25

applications, for instance: genetics (Jakobsson and Rosenberg , 2007; Altuvia

26

et al. , 2005), medicine (Khanmohammadi et al. , 2017), marketing (Ngai et

27

AC

Unsupervised learning techniques have paramount importance when unla-

∗ Corresponding author. Preprint to [email protected] Expert Systems with Applications May 30, 2017 Emailsubmitted addresses: (Ariel E. Bay´ a), [email protected] (M´ onica G. Larese), [email protected] (Rafael Nam´ıas)

ACCEPTED MANUSCRIPT

al. , 2009), remote sensing (Tarabalka et al. , 2009), web mining (Srivastava

2

et al. , 2000), text mining (Dhillon and Modha , 2001; Tseng et al. , 2007)

3

and image segmentation (Shi and Malik , 2000; Mur et al. , 2016), to name a

4

few. In fact, image segmentation and unsupervised learning share a bond, as

5

described in (Shi and Malik , 2000; Arbelaez et al. , 2011; Comaniciu and Meer

6

, 2002; Felzenszwalb and Huttenlocher , 2004). Consequently, common issues

7

arise. On the one hand, from the unsupervised learning viewpoint, how many

8

clusters can be found in a set of patterns. On the other hand, from the image

9

processing perspective this topic could be formulated as: how many objects can

10

be identified in an image. The perspective from each community thoroughly

11

diverges. On the one side, papers from the image processing community (Alata

12

and Ramananjarasoa , 2005; Bergasa et al. , 2000; Chen et al. , 2005) rely en-

13

tirely on image processing methods to deliver the final result while, on the other

14

side, publications with a machine learning approach (Herbin , 2001; Stanford

15

and Raftery , 2002; Aue and Lee , 2011) make focus in validation processes.

AN US

CR IP T

1

Clustering validation research is usually applied to bioinformatics or biolog-

17

ical problems (Handl et al. , 2005; Yeung et al. , 2001; Bay´ a and Granitto ,

18

2013; Rousseeuw , 1987). Mostly because biological data is high dimensional

19

and has few samples compared to its dimensionality, data clustering and cluster-

20

ing post-processing are essential to analyze this kind of data. On the contrary,

21

image processing is a much more mature discipline and it has more techniques

22

available for post-processing. Nevertheless, most of them require handcrafted

23

solutions that cannot be directly applied to different domains, whereas solutions

CE

PT

ED

M

16

24

based on learning methods try to keep customization to a minimum to enable

25

portability across knowledge domains.

AC

26

We consider image segmentation as a special type of unsupervised learn-

27

ing problem. Hence, our solution combines clustering and validation to achieve

28

segmentation/clustering without using any further image post-processing. To

29

achieve our goal, we modify a clustering validation method, Clustering Stabil-

30

ity (CS) (Lange et al. , 2004), to partition/segment: textures, natural images

31

and simulated 3D MRI data in an automated fashion. First, we extract color 3

ACCEPTED MANUSCRIPT

1

data driven method, i.e., a non-perceptual procedure, able to find the optimum

2

number of segments within an image according to a new color-texture repre-

3

sentation. Our strategy merges a texture descriptor (Ojala , 2002) with color

4

CR IP T

and texture features to create a new image representation. We introduce a

5

it usually requires to adjust several parameters to find a proper segmentation.

6

Although CS provides an unsupervised score that automatically allows to find

7

the number of clusters/segments, we extended this score to also evaluate a set

8

of representations coming from diverse parameters settings. Therefore, we our

9

method can simultaneously perform image segmentation and select the most

10

appropriate data representation as well.

11

AN US

information. Finding an useful color-texture representation is demanding, since

12

ysis in the image processing literature. Ilea and Whelan (2011) present a

13

comprehensive review including most relevant methods up to 2011. The appli-

14

cations considered include: object recognition, scene understanding, image in-

15

dexing and retrieval, among others. Up to our knowledge, current unsupervised

16

segmentation methods rely heavily on ellaborated or taylored image processing

17

operations, even those based on pattern recognition algorithms. For example,

18

ED

M

There is a wide range of segmentation examples based on color-texture anal-

19

dominant colors and merge them with texture features obtained with Gabor

20

filters. The aim of the proposed features is to detect homogeneous regions in

21

the image and integrate them by the use of an adaptive clustering algorithm

22

enforcing the spatial continuity by minimizing an energy function. Kim and

23

CE

PT

Ilea and Whelan (2008) proposed the use of RGB and YIQ color spaces to find

24

with a clustering strategy aimed for texture feature extraction. Finally, the

25

authors merged texture features and RGB color information to segment the

26

images with a graph cut algorithm that minimized an energy function. This

27

scheme involving color-texture features extraction and clustering has also been

28

successfully used by Yang et al. (2013) and de Luis Garc´ıa (2008). This kind

29

of framework involves a new type of color-texture features, which transforms

30

the original image into a new representation used as a proxy. Then, a cluster-

31

AC

Hong (2009) explored the use of Gabor filter banks as well, by combining it

4

ACCEPTED MANUSCRIPT

ing algorithm can be applied to divide the new representation leading to the

2

segmentation of the original image. A negative aspect of using color-texture

3

features is the existence of free parameters, which all were set by hand. A more

4

desirable scenario would be to have an unsupervised score allowing to compare

5

multiple solutions arising from diverse settings in order to select the best one.

6

Our work not only modifies CS to select the number of clusters/segments in an

7

image, but also demonstrates that CS can be used to detect the optimal settings

8

required for a pre-processing stage like color-texture feature extraction.

CR IP T

1

More recently, Mur et al. (2016) presents a method to determine the opti-

10

mal number of clusters using spectral clustering, the Silhouette validity index

11

and local scaling for image segmentation purposes. They test their method

12

on synthetic data and gray scale images. The computational cost required for

13

methods such as spectral clustering is high in both time and memory due to re-

14

quirements necessary for the distance matrix and eigenvectors computation and

15

storage. The authors propose low computational cost versions of their method

16

in order to achieve more practical efficiency.

M

AN US

9

In brief, we not only introduce some major modifications to CS, but also

18

discuss a new representation for the color and texture features in order to achieve

19

memory efficiency. In our experiments on texture and color natural images,

20

we are able to automatically find the best image segmentation and texture

21

parameters simultaneously. We also test and discuss our method when applied

22

to synthetic 3D MRI data. Additionally, our method provides a stability score

23

profile which allows to easily inspect alternative segmentation solutions.

CE

PT

ED

17

The paper has four sections. After the introduction, in Section 2, we describe

25

the color and texture feature extraction procedure and how we adapt CS for

26

image segmentation purposes. Then, Section 3 presents experimental results

AC

24

27

on texture, real images and MRI data. An additional file is attached where we

28

test an alternative texture descriptor and provide more experimental results.

29

Finally, Section 4 presents our conclusions.

5

ACCEPTED MANUSCRIPT

2. Image segmentation by means of CS

1

CS can be used to divide an image according to its information content. In

2

this paper, we use CS to process gray scale intensities, color and texture images

3

from an image as a

4

describes any pixel information content.

5

For instance, for RGB color space pi,j ∈ {R, G, B}. However, the method k

6

is not restricted to any particular color space. Moreover, a pixel p

could be

7

described by any number of features including texture information as well. Thus,

8

pi,j = {color features} + {texture features}, where {color features} describes

9

the color from pixel (i, j), {texture features} describes texture information from

10

pixel (i, j) and the + symbol concatenates both feature vectors into a new vector

11

representation.

12

i,j

vector, thus, p

=

i,j i,j (pi,j 1 , p2 , ...., pn )

CR IP T

both separately and combined. We consider a pixel p

i,j

AN US

i,j

i,j

13

task is to automatically divide them into groups. The pixel information content

14

is based on gray scale intensities, color, texture, or both color and texture. To

15

M

We consider an image I as a collection of pixels p , and our unsupervised

16

tering plus validation. First, block processing is applied to the image, dividing

17

it into non-intersecting patches from which features are extracted. The second

18

step, clustering and validation, is usually treated as two separate procedures

19

(Handl et al. , 2005), but since CS performs both simultaneously we prefer to

20

PT

ED

cope with this issue we divide it in two steps: 1) feature extraction and 2) clus-

21

following subsections.

22

CE

consider them as a single operation. We explain steps 1 and 2 in detail in the

2.1. Feature extraction

23

24

tures. In our experiments, we use the following features: gray scale intensities,

25

color features, texture features and a combination of color and texture. We use

26

RGB to represent color images. We pre-process them as follows: i) we remove

27

each channel (R, G, B) mean and ii) we apply PCA whitening to step (i) out-

28

put. Whitening will standardize deviation along principal components to 1. We

29

AC

CS is a versatile method which can be used to process different kinds of fea-

6

ACCEPTED MANUSCRIPT

standardize data in order to combine color and texture features, which can have

2

different scales. Thus, we avoid one set of features dominating the other. Also,

3

since most of our experiments compare color and color + texture we use the

4

same color normalization in both cases to be consistent. The following section

5

2.1.1 explains the texture feature extraction procedure as well as discussing deep

6

into details the data normalization.

7

2.1.1. Texture features

CR IP T

1

For texture coding we selected a descriptor based on a 2D histogram orig-

9

inally developed by Ojala (2002). The authors use the discrete distribution LBP RI V AR ,

AN US

8

where LBP RI is a gray scale rotation invariant (RI) local bi-

10

HRI =

11

nary pattern (LBP) and V AR measures the contrast of local image texture,

12

which is invariant to rotation but not to scale. Locality is addressed in terms

13

of pixel neighborhoods, which, in our paper, we set using windows of size

14

w ∈ {11, 15, 21, 25, 31, 35, 41, 45, 51}.

However, HRI descriptive power is not enough for our segmentation purpose

16

because LBP RI has few states, which are targeted for specific features, i.e.

17

line, corners, etc. Thus, we decided to change LBP RI to LBP N RI , a non-

18

rotation invariant (NRI) LBP. NRI-LBP provides a higher degree of texture

19

orientations, yet, this improvement comes at a high price. While an RI-LBP

20

with an 8-point neighborhood and radius 1 requires of 9 codes to describe a

21

texture, a NRI-LBP needs 256 codes. In spite of this difference the method to

22

find this descriptor is analogous to the one described by Ojala et al. to model

PT

ED

M

15

23

LBP RI V AR .

24

of size 256 × 9 by using HN RI =

25

(M × N × 256 × 9) to describe the complete texture for an image of size M × N .

LBP N RI V AR .

To do so, we use a tensor of size

AC

CE

As a result, we can describe the texture at pixel (i, j) as a 2D histogram

26

This tensor representation can potentially consume a great amount of memory.

27

For example, a case where M × N × 1 Byte = 500 kilobytes would need over

28

1 terabyte of memory to represent the complete texture tensor. Therefore, to

29

reduce memory requirements we used patches of size L × L rather than pixels.

30

In this work we used L = 5, but higher values can be used instead depending 7

CR IP T

ACCEPTED MANUSCRIPT

Figure 1: Example of image patches for texture representation.

1

size.

2

AN US

on the available memory resources. Therefore, henceforwards, L = 5 for patch

3

image in L × L non-overlapping patches. Then, within each patch, we compute

4

the L × L × 256 × 9 NRI-LBP descriptor. Roughly, using patches we reduce

5

the memory requirement by 3 orders of magnitude. Figure 1 shows an example

6

using a 2 × 2 pixel (L = 2) patch. The notation Iu,v denotes that pixel I

7

belongs to patch u and is the pixel number v from the patch. The code using

8

M

To model an image texture using the HN RI LBP descriptor we divide the

10

is the number of pixels.

12

u,v HN RI

ED

(u, v) is equivalent to (i, j), yet more appropriate in this context. The NRIP u,v 1 u u LBP histogram HN RI for patch u is calculated as: HN RI = S v HN RI , where

is the histogram corresponding to pixel v within patch u and S = L × L

9

11

13

the texture. Afterwards, we require to provide CS with the complete texture

14

information, even if the information is the same, i.e., all L × L pixels v within

15

patch u have the same

Thus, we need to compress the texture informa-

16

tion to keep the most relevant part of it. Therefore, we post-process the texture

17

u descriptors (HN RI ) with PCA and we only save a small number of components.

18

u HN RI .

AC

CE

PT

The use of patches helps to reduce memory usage, but only while modeling

This process requires a few steps: i) vectorize histograms from all groups to form   ×N a (M L×L ) × (256 × 9) matrix and ii) apply PCA to them. The final texture

19

20

descriptor is a vector of P components, where P  (256 × 9). Finally, we can

21

describe each pixel (i, j) texture by a P component vector, which only requires

22

to map the coordinates (u, v), coding the number of patch, to the pixel number

23

8

ACCEPTED MANUSCRIPT

1

(i, j) according to the image coordinate system. In any case, this transformation

2

is trivial. To describe both texture and color, first we concatenate color and texture

3

5 6

descriptors, then we normalize the data columns to zero mean and unit standard q 1 deviation, and at last, we scale the P texture components by P , thus, P

texture features weight as a single color feature/channel.

CR IP T

4

7

Finally, we consider an additional texture descriptor based on Gabor filters

8

in order to verify that CS and the NRI-LBP texture descriptor are not coupled

9

and do not give a biased result. The Gabor texture descriptor implementation and derived experimental results are included in an additional file.

11

2.1.2. Sampling

AN US

10

The sampling procedure divides the image into two disjoint sets, A and

13

B, necessary for CS (see Section 2.2.1). Figure 2 shows a simplified diagram

14

illustrating this process. The image is divided into non-overlapping L×L patches

15

and features previously computed are extracted for random pixels from the

16

blocks and assigned to one of the two sets, A or B. By these means, we preserve

17

the original pixel spatial distribution in both sets. We use pi,j to refer to the

18

pixel p at position (i, j), where i is a row and j is a column. Pixel pi,j has

19

probability P r = 0.45 to be added to set A or B respectively. The sampling

20

step starts by generating a random value pr0 , where 0 ≤ pr0 ≤ 1. Then, if

PT

ED

M

12

21

the proposition P r < pr0 is true, we add pi,j to Set A. On the other hand, if

22

P r ≤ pr0 < 2×P r is true we add pi,j to Set B. We drop the pixel if pr0 ≥ 2×P r. We repeat this procedure for each pixel pi,j+1 to the last column. Then we

24

proceed in the same fashion with the next row, i.e. pi+1,j . We apply the

25

sampling to every image patch as we show in Figure 2.

AC

CE 23

26

CS repeats this sampling step several times, and A and B sets cardinality

27

does not change. However, pixels selected by this method ensure that new

28

instances of sets A and B consider pixels with different spatial coordinates. But

29

at the same time, both sets keep local information since pixels are extracted

30

from small local patches. This downsampling method prevents the inclusion of 9

AN US

CR IP T

ACCEPTED MANUSCRIPT

Figure 2: The diagram shows the sampling procedure for feature extraction using (L = 5) and (P r = 0.45) as explained in Subsection 2.1.2. A fraction (P r) is assigned to each set (A and B) and the remaining points (10%) are dropped. In this example, the pixels in green are

M

assigned to Set A, the pixels in white are assigned to Set B and the pixels with a red cross are ignored.

1

information. Instead, we force both sets to have a fraction of P r = 0.45 of the

2

pixels per patch.

3

PT

ED

long pixels strings into either set, which could deprive the other one from local

2.2. Clustering and validation

4

5

idation framework which is able to find the best data partition. An advantage

6

from this validation algorithm is that it can use any clustering technique. Apart

7

from this, its computational complexity increases with the number of iterations,

8

and its spatial complexity is linear. However, due to the intrinsic data size, some

9

applications need both fast clustering and validation algorithms. For instance,

10

MRI images. Hence, we have only considered a fast clustering algorithm like k-

11

means within our modified version of CS. Figure 3 presents a simplified diagram

12

AC

CE

CS uses a clustering algorithm to divide the image into clusters within a val-

10

ACCEPTED MANUSCRIPT

of the clustering and validation stage. The clustering algorithm is actually re-

2

lated to other elements within the validation method. Therefore, the validation

3

has to be performed according to a set of rules and considerations which are

4

later described in the next section. Finally, since clustering and validation are

5

related we consider them as a single operation.

CR IP T

1

Figure 3 summarizes the pipeline behind the clustering and the number

7

of clusters selection. First, the data needs to be partitioned into two sets A

8

and B. While the original method uses a random sampling dividing data in

9

halves, we implemented the sampling strategy from Subsection 2.1.2. Then,

10

we clustered each set and trained a classifier using the labels found by the

11

clustering algorithm resulting in two sets of labels: i) clustering labels for Set A

12

and predicted labels for Set B, and ii) clustering labels for Set B and predicted

13

labels for Set A. Finally, the labels resulting from clustering Set A are compared

14

with the predicted labels of Set A, which are found by using a classifier trained

15

with Set B. Analogously, we perform the procedure for Set B. The outcome

16

from this comparison is a score in the range [0, 1], where 0 stands for perfect

17

match and 1 for non match. This process is repeated a number of times using

18

different samples, i.e., applying the sampling strategy from subsection 2.1.2 R

19

times. The final outcome from CS is the mean score considering R iterations.

20

Fellowing, we describe CS in more detail. The reader may be referred to Lange

21

et al. (2004) for further information on this procedure.

22

2.2.1. Clustering stability

PT

ED

M

AN US

6

CS starts by dividing the original data in two sets, A and B. A major consid-

24

eration is that each set needs to replicate the original data distribution to ensure

25

the CS success. When this condition is met, we can use both sets to infer the

AC

CE 23

26

original data clustering structure. We use the procedure from Subsection 2.1.2

27

as sampling procedure, which saves spatial relations and pixel distribution.

28

CS requires to couple a clustering algorithm and a classifier. The clustering

29

algorithm, for instance, k-means (Macqueen , 1967), is used to explore and to

30

find groups in the data. Thus, the number of clusters used to run k-means 11

CR IP T

ACCEPTED MANUSCRIPT

Figure 3: Pipeline of clustering stability procedure. The figure describes the complete validation process starting at Feature Extraction described in Section 2.1. The output from the

feature extraction process is sampled according to the procedure explained in Section 2.1.2.

the validation score named STI index.

AN US

Next, the clustering stability method from Section 2.2.1 is applied, which finally calculates

1

{L2 , L3 , ...., Lkmax }, where Lj are the labels found by k-means when the number

2

of clusters is set to j. The labels Lj are then used to train (kmax − 1) classi-

3

fiers, one per clustering partition. A good practice is to select the most simple

4

classifier able to learn the clustering partition. However, selecting a classifier

5

not able to learn the clustering partition will drive CS to a bad solution. In

6

case of k-means, any linear classifier will yield a good solution. Yet, to find

7

ED

M

has to be varied from 2 to kmax . As a result, there are (kmax − 1) partitions

8

dates simultaneously with a fixed clustering algorithm. We selected the Linear

9

Discriminant Analysis Classifier (LDAC) (Hastie et al. , 2003). We predict the

10

labels using LDAC from the set which was not used for clustering. Thus, if Set

11

A is clustered then its labels are used to train (kmax − 1) classifiers. Finally,

12

Set B labels are predicted with each classifier. CS uses both sets for clustering

13

and prediction, resulting in two sets of labels for A and B, one from clustering

14

and the other one from prediction. Figure 3 shows a flow diagram of the CS

15

AC

CE

PT

which classifier works best, according to CS, we would need to test the candi-

procedure.

16

We name the sets of labels resulting from clustering Set A and B as La and

17

Lb , respectively. The elements within L∗ are L∗ = {L(2,∗) , L(3,∗) , ...., L(kmax ,∗) },

18

where ∗ is either a or b. We use the clustering label sets to train two groups of

19

classifiers: i) the first classifier group uses La and Set A and ii) the second uses

20

12

ACCEPTED MANUSCRIPT

Lb and Set B. Thus, each group trains a set of (kmax −1) classifiers C(i,∗) , where

2

i = 2, 3, ..., kmax indicates the number of clusters in the partition and ∗ = {a, b}

3

indicates the set. For example, C(2,a) is a classifier that learned Sets A two

4

clusters partition.

CR IP T

1

After the clustering/training stage, we have a set of clustering labels and

6

classifiers {La , Ca } for Set A and one for Set B {Lb , Cb }. CS uses the information

7

from these two sets to derive a score that can infer the number of clusters existing

8

in the original data. The score used by CS compares the similarity between two

9

partitions through a classifier. This concept is known as transference. Next,

10

we provide a transference example using Set A as origin point. Let us consider

11

two partitions of Set A and B with (k = 2), i.e., L(2,a) and L(2,b) . First, we

12

train a classifier using Sets B labels L(2,b) and name it C(2,b) . Next, we use

13

C(2,b) to predict Sets A labels, which results in C(2,b) set that partitions Set

14

A according to the structure learned from Set B. Finally, we need to measure

15

the similarity between C(2,b) and L(2,a) . Both sets partition Set A, however, we

16

ignore the correspondence among the labels from L(2,a) and C(2,b) . Therefore,

17

we need to solve a matching problem so we can measure this relationship. The

18

correspondence can be found by applying the Hungarian algorithm (Kuhn ,

19

1955). Once the sets are matched, we can calculate the similarity between sets.

20

Figure 3 refers to this value as Score A and Score B. These scores measure

22

(2,a)

ED

M

(2,a)

the adjusted error rate (Ek,a ). We first calculate the error rate by: Ek,a = P 1 (k,a,i) , where Ma = #(Set A), L(k,a,i) is point i clustering i {1}{L Ma 6=C } (k,a,i)

(k,a,i)

label and C(k,b)

CE

23

(2,a)

PT

21

AN US

5

(k,b)

is point i predicted label. The score Ek,a is the error using Set

24

A as origin point, hence, we could also find an analogous score Ek,b using Set

25

B as origin point. The adjusted error rate formula is Ek,a =

Ek,a −Ek,a , max −E Ek,a k,a

where

27

value of Ek,a under the null hypothesis. For further information, there is an

28

interesting discussion on normalization and correction for chance by Vinh et al.

29

(2009).

AC

26

max Ek,a

30

is the maximum numerical value of the index and Ek,a is the expected

We use the notation Ek to refer to the Stability Index (STI). We calculate

13

ACCEPTED MANUSCRIPT

the STI index by repeating the clustering/training stage R times with Set A as origin point. Thus, we end up with R error measures

z (Ek,a ),

1

where k =

2

{2, 3, ..., kmax } refers to the number of clusters, and z = {1, 2, ..., R}, to the

3

4

CR IP T

z number of repetitions. Likewise, there is an error measure (Ek,b ) arising from

6

value of STI index (Ek ) for each k = {2, 3, ..., Kmax } and the lowest STI (Ek )

7

points to the “best” solution. In this context, the “best” solution is the one

8

having a consistent low error, since outliers tend to be overlooked. Lower errors

9

indicate that for a given (k = c∗), samples of the original data can model each

10

clustering structure.

11

AN US

the use of Set B as origin point, which we can use to boost the statistics. Finally, P z 1 z CS uses the (2 × R) values to find: Ek = 2×R z Ek,a + Ek,b , where there is a

3. Experimental results

5

12

13

segmentation, ii) color and texture segmentation of natural images and iii) syn-

14

M

In this section we analyze the obtained results for three tasks: i) texture

15

parts are composed by textures and natural images without ground truth, thus

16

we can only evaluate (perceptually) the results from CS. The images used in

17

part i) and ii) are of public domain and they can be requested by e-mail to

18

the authors. For these experiments we also use The Berkeley Segmentation

19

PT

ED

thetic 3D brain MRI segmentation. The experimental results from the first two

20

we used synthetic 3D MRI images from the Simulated Brain Database1 created

21

with BrainWeb (Cocosco et al. , 1997). For these images, the true labels are

22

available, allowing us to compare them with the clustering labels found by CS

23

to quantitatively evaluate the segmentation performance.Table 1 describes the

24

main features of all the data used in this paper. For 2D images, the number of

25

patterns is the total number of pixels.

26

AC

CE

Dataset and Benchmark (Martin , 2001). On the other hand, for the third task

On the other hand, for MRI data the majority of voxels are zero (back1 http://www.bic.mni.mcgill.ca/brainweb/

14

27

ACCEPTED MANUSCRIPT

ground). Since the task in this case consists in segmenting gray matter, white

2

matter, cerebrospinal fluid, edema and tumor, the background voxels are dis-

3

carded, highly reducing the number of patterns to be processed. To segment the

4

images, we apply CS with k-means as clustering algorithm and LDAC as classi-

5

fier, obtaining scores ranking the clustering partitions that we use to select the

6

best partition. Moreover, we use CS not only to select a clustering/segmentation

7

solution, but also to set the window size required for HN RI to calculate the tex-

8

ture features. Thus, we consider several w × w windows sizes without any a

9

priori knowledge about the best window size. In our experiments, “the best

10

window size” is the one that minimizes partitions scores. Finally, to build the

11

texture feature vector we set P = 12 in all experiments. The first two tasks,

12

texture segmentation and segmentation of natural images, simultaneously select

13

the best partition/segmentation solution and the window size. For this purpose,

14

we compare STI vs. Number of clusters plots. The legend of these plots are

15

the conditions we are testing, these are the window sizes used to calculate the

16

NRI-LBP texture descriptor and the absence of texture, which we name “color”.

17

The idea is to find the “best” window size according to the STI index. However,

18

it could be possible that the use of texture worsens the segmentation. Hence, we

19

test color alone too. CS selects solutions near the absolute STI minimum using

20

a threshold (T=0.05). We call these solutions candidates. Candidate solution

21

tells us which conditions, i.e., window size and/or absence of texture, have a

22

minimum near the global minimum and where it happens, which amounts to

23

the value of k. We use an orange circle to point out the candidate solutions,

CE

PT

ED

M

AN US

CR IP T

1

the remaining parameters to be set are the number of times CS repeats the

25

clustering / learning cycle (R) and the maximum number of clusters (Kmax )

26

CS looks for. For all experiments we set R=30 and Kmax = 10.

AC

24

27

The final task, synthetic 3D T1 MRI data segmentation, compares the per-

28

formance of different clustering validation methods for finding the best partition.

29

Since the type of tissue simulated in the MRIs is known, voxels belong to one of

30

the following classes: gray matter, white matter, cerebrospinal fluid, edema and

31

tumor. The experiment consists in applying CS, the WB (Zhao et al. , 2009) 15

ACCEPTED MANUSCRIPT

Table 1: Main features of the considered datasets. (*) Foreground voxels considered for segmentation. Image

Type

Size (M × N )

Synthetic texture

Gray scale

256 × 256

Wood texture

RGB

296 × 950

281200

Wool texture

RGB

863 × 1383

1193529

Peppers

RGB

757 × 1051

795607

Limes

RGB

1200 × 1600

1920000

Landscape A

RGB

400 × 600

240000

Landscape B

RGB

1030 × 1800

1854000

BR 01

T1

256 × 256 × 181

2021427 (*)

BR 02

T1

256 × 256 × 181

2062629 (*)

BR 03

T1

256 × 256 × 181

2062865 (*)

BR 04

T1

256 × 256 × 181

2016961 (*)

BR 05

T1

256 × 256 × 181

2020970 (*)

BR 06

T1

256 × 256 × 181

2062811 (*)

No. of patterns (pixels/voxels)

AN US

CR IP T

65536

1

then using MRI true labels to rate the segmentation achieved by validation

2

methods. Partition assigns each voxel to a cluster, hence, cluster membership

3

can be contrasted against true labels. However, before this can be done, cluster-

4

ing labels and true classes need to be matched using the Hungarian algorithm

5

(Kuhn , 1955). We assess the quality of the method by computing the accuracy

6

ED

M

and CH (Calinski and Harabasz , 1974) indexes to find the best partitions and

7

match between clusters and classes and 0 means that no matches were made.

8

As part of the results, we also report the number of clusters detected by each

9

PT

of the segmentation. This figure lays within [0, 1], where 1 stands for perfect

validation method.

10

11

In this work we consider texture features for texture images, color and texture

12

features for natural images, and gray scale intensities for MRI data.

13

3.1. Texture experiments

14

AC

CE

We show how CS can be used in combination with different types of features.

We used the three images from Figure 4 as texture segmentation cases of

15

study. The test involved applying CS with k-means+LDAC to find the image

16

partitions, which we present as a gray scale mask. We also used CS as the tune

17

up method to select a parameter for texture extraction. Following Section 2.1

18

16

(a)

(c)

AN US

(b)

CR IP T

ACCEPTED MANUSCRIPT

Figure 4: Texture images. (a) Wood. (b) Synthetic. (c) Wool.

1

we had several window size candidates w = {11, 15, 21, 25, 31, 35, 41, 45, 51} to

2

extract texture features by means of HN RI . Thus, our objective was to find the

3

best image partition and window size simultaneously.

Figure 5 panels present three typical CS profiles, which in this case provide

5

a solution to segment the texture images. The (scores) profile x-axis reports

6

the number of clusters and its y-axis reports the score values (STI) ranking the

7

clustering partitions. The panels include several STI evolutions, each of them

8

illustrates the clustering scores for a given window size. For the sake of clarity,

9

we include only the most relevant windows sizes, i.e. the ones closest to the

10

global minimum. Also, we use an orange circle to point out relevant solutions.

PT

ED

M

4

In each case the best solution corresponds to the minimum STI: w = 35 and

CE

11

k = 2 for wood texture, w = 21 and k = 5 for synthetic texture and w = 21 and

13

k = 3 for wool texture. Figure 6 shows the gray scale masks obtained for the

14

three textures, and an extra wood texture mask we consider relevant to discuss.

AC

12

15

In the wood texture example, a window size of w = 35 leads to the best

16

partition. The mask corresponding to this partition is shown in Figure 6(a).

17

This mask divides the image in two clusters according to the orientation of their

18

elements: vertical or horizontal. The extra wood texture mask corresponds to

19

k = 4 and w = 25 in Figure 5(a) (also marked with an orange circle). This 17

AN US

(a)

CR IP T

ACCEPTED MANUSCRIPT

(b)

(c)

Figure 5: Scores obtained for the texture images and different window sizes. Orange circles

M

point out relevant solutions. (a) Wood. (b) Synthetic. (c) Wool.

1

orientation. Figure 6(b) has four clusters one per element: planks and spaces

2

ED

partition divides the image in planks and spaces between planks according to

times two orientations, horizontal and vertical. The mask shows that some parts

3

of the planks are mislabeled as if they were spaces.

4

PT

The following mask, shown in Figure 6(c), divides the synthetic texture in

5

6

though the error tends to increase near the borders due to the coarseness of the

7

fabric.

8

CE

five pieces, one corresponding to each texture. The partition is quite accurate,

9

pattern difficulty arises from the knitting trace pseudo regularity. To rate the

10

AC

Finally, the mask for the wool texture is displayed in Figure 6(d) where the

quality of the mask, we analyze three sub-patterns: the texture has two levels,

11

the white cluster from Figure 6(d) being the shallow section of the texture.

12

The rest is formed by two sub-patterns with variable height crossing each other.

13

These sub-patterns are shown as two clusters in Figure 6(d): one black and one

14

gray. The black cluster corresponds to the pattern crossing from above whereas

15

18

CR IP T

ACCEPTED MANUSCRIPT

(a)

AN US

(b)

(c)

(d)

M

Figure 6: CS best segmentation masks for the texture images. (a) Wood image with w = 35 and k = 2. (b) Wood image with w = 25 and k = 4. (c) Synthetic image with w = 21 and

ED

k = 5. (d) Wool image with w = 21 and k = 3.

the gray cluster corresponds to the pattern crossing from below. Figure 6(d)

2

also shows that areas near sub-patterns crossings have higher error.

3

3.2. Color and texture experiments

PT

1

As in the previous subsection we use CS with k-means for clustering and

CE

4

LDAC for learning. We pursue the segmentation of natural images, but now

6

considering color and texture features together. We include texture information

7

as a way to improve cluster detection. The color and texture features merge

AC

5

8

requires to tune the window size w for HN RI . As previously described, we select

9

the best partition and w simultaneously. Moreover, we include a supplemen-

10

tary condition for comparison purposes, which is the usage of color alone. We

11

use this setting to check if our combined set of features actually improves color

12

segmentation. Our plots present the STI values for the partitioning with dif19

ACCEPTED MANUSCRIPT

1

reports the number of clusters and the y-axis reports the STI values ranking

2

the clustering partitions. The natural image segmentation test aims at selecting

3

simultaneously the best clustering partition and window size, or no window in

4

CR IP T

ferent number of clusters. We draw score profiles as before, where the x-axis

case that color alone is the best choice for segmentation. These experiments use

5

the following window sizes w = {31, 41, 51, 61, color}, where numbers refer to

6

sizes in pixels and color means color segmentation.

7

8

a white background. Figure 7(a) and Figure 13(a) show an image of the pep-

9

pers and the corresponding STI profiles considering various window sizes, re-

10

spectively. In this case the best partition has 4 clusters, independently of the

11

window size. Figures 7(b) and (c) show the masks resulting from the usage of

12

color features alone, and texture and color combination, respectively. It can be

13

noticed that there are subtle differences between the two masks. Both masks

14

divide the image in background (which includes the peppers shadow) plus one

15

M

AN US

Peppers: We consider a simple image of three peppers casting a shadow over

16

pepper also holds the green stalks from each pepper. In this example, all five

17

conditions for w produce the same solution, as it can be depicted in Figure

18

13(a).

19

ED

cluster per pepper (red, yellow and green). The cluster including the green

20

captured against an uneven background. Figure 13(b) shows the STI profiles

21

ranking partitions and conditions. Most features have the same STI minimum

22

CE

PT

Limes: Figure 8(a) consists of a photo of limes plus some slices of the same fruit

23

higher STI value. There are also local minimums at (k = 7, color) and (k = 8,

24

w = 41). Figures 8(b) and (c) show the masks dividing the image into 4 and 7

25

AC

value pointing at 4 clusters with exception of w = 41, which has a slightly

clusters respectively, using color information only, and Figure 8(d) displays a 4

26

clusters mask resulting from texture and color segmentation. The two 4 clusters

27

masks, from Figures 8(b) and (d), differ only in small details. For example,

28

in Figure 8(d) we can see the lime pulp as a single cluster, but Figure 8(b)

29

considers some white pulp bits as part of the background. The white pulp bits

30

20

(a)

(c)

AN US

(b)

CR IP T

ACCEPTED MANUSCRIPT

Figure 7: Peppers. (a) Original image. (b) Color features. (c) Color+texture features.

and the background have different textures, therefore, texture features increase

2

the distance between them. As a result, the lime pulp is included in a single

3

cluster. However, there are non-related parts of the image that are still grouped

4

together. For example, part of the white membrane between the pulp and skin

5

is clustered with the background in Figures 8(b) and (d), and in less amount in

6

Figure 8(c). Unlike the white pulp bits, this membrane is a wide white strip with

7

uniform texture, which is very similar in color and texture to the background.

ED

M

1

Finally, Figure 9 presents a detailed view of the partitions from the mask

9

previously shown in Figure 8(c) corresponding to the local minimum at (k = 7,

10

color). This solution is an over-fragmented version of the 4 clusters parti-

CE

PT

8

tion from Figure 8(b). We can reconstruct a similar solution to the one at

12

(k = 4, color) by merging the clusters (2 ⊕ 3), (4 ⊕ 5) and (6 ⊕ 7). Thus,

13

{1, {2, 3}, {4, 5}, {6, 7}} results in the partition at (k = 4, color). In spite of

AC

11

14

the fact that solutions at (k = 7, color) and (k = 4, color) are related, STI

15

values are not similar. A higher STI tells us that the divisions between clus-

16

ters are not clear. To understand this, consider a dataset without clear clusters

17

boundaries, which we sample and cluster repeatedly. If clusters boundaries are

18

not clear, patterns near the borders can change membership from one iteration 21

CR IP T

ACCEPTED MANUSCRIPT

(b)

AN US

(a)

(c)

M

(d)

Figure 8: Limes. (a) Original image. (b) k = 4, color features. (c) k = 7, color features. (d)

ED

k = 4, color+texture features.

1

iterative clustering+learning process from Figure 3 to a dataset without clear

2

cluster boundaries. In this case, we find high STI scores. On the contrary,

3

if we consider moving border points closer toward its centroid by compressing

4

CE

PT

to the next. Following this line of thought, let’s consider that we apply the

5

contributes to improve cluster division, which makes partitions easier to learn.

6

AC

the cluster uniformly, this action will reduce the previous STI values because it

Landscape A: The outdoor scene from Figure 10(a) has three distinctive ele-

7

ments: grass, sky and a signpost. The resulting STI profile is displayed in

8

Figure 13(c). We select solutions with k = {2, 4}, which have almost the same

9

STI, plus a third solution with k = 6 and higher STI. Figures 10(b-e) show

10

masks dividing the image into 4, 6, 2 and 4 clusters respectively. Panels (b)

11

22

AN US

CR IP T

ACCEPTED MANUSCRIPT

Figure 9: Detailed view of the clusters from Figure 8(c) (k = 7, color features).

and (c) are the result of CS and color features alone while panels (d) and (e)

2

result from CS using color in combination with texture. We can see that both

3

4 clusters solutions, Figures 10(b) and (e), are quite different. For example, the

4

lines dividing the sky are not similar and the signpost deer is mostly included

5

in the same cluster as the dark part of the sky in Figure 10(b).

M

1

As we discussed above, over-fragmented solutions without clear cluster di-

7

vision are penalized by CS. Figure 10(c) is an example of this, the sky and the

8

grass are divided at the STI expenses. Finally, Figure 10(d) presents a 2 clusters

9

solution dividing the sky from the rest of the picture. This solution has a low

10

STI value similar to the one of Figure 10(e), thus we have two possible solutions

11

with minimum STI values.

CE

PT

ED

6

Landscape B: Figure 11(a) displays a coastal scene with three major compo-

13

nents: sea, terrain and vegetation. There is also a smaller component that

AC

12

14

could form a fourth cluster: a white foam layer over the water. STI profiles

15

shown in Figure 13(d) point to solutions with k = 3 clusters. We also analyze

16

a local minimum at (k = 7, color) with a higher score. The masks partitioning

17

the image are displayed in Figures 11(b-d), where panels (b) and (c) correspond

18

to 3 and 7 cluster partitions using color features and panel (d) to a 3 cluster

23

AN US

CR IP T

ACCEPTED MANUSCRIPT

M

(a)

PT

ED

(b)

(d)

(c)

(e)

CE

Figure 10: Landscape A. (a) Original image. (b) k = 4, color features. (c) k = 6, color

AC

features. (d) k = 2, color+texture features. (e) k = 4, color+texture features.

24

(a)

(c)

AN US

(b)

CR IP T

ACCEPTED MANUSCRIPT

(d)

Figure 11: Landscape B. (a) Original image. (b) k = 3, color features. (c) k = 7, color

1

M

features. (d) k = 3, color+texture features.

partition using color+texture features.

The solution with higher STI score, (k = 7, color), over-fragments the im-

3

age: terrain, sea and vegetation are divided each in two clusters. The ocean

4

foam is assigned to a separate cluster, which includes also a small coastal line

5

portion made of light gray stones. On the contrary, solutions with 3 clus-

6

ters divide the image into: sea, terrain and vegetation. Moreover, the seg-

7

mentation based on color+texture features achieves identical results for w =

CE

PT

ED

2

{31, 41, 51, 61}. However, there are notable differences between color (Figure

9

11(b)) and color+texture based segmentation (Figure 11(d)). We include Fig-

10

ure 12 to make easier the comparison of both masks. Figure 12(a) shows three

11

panels, where each of them corresponds to a cluster: sea, terrain and vegetation,

12

that were obtained using color features. Analogously, Figure 12(b) displays the

13

same three clusters resulting from color+texture features. A comparison of both

14

left panels, the sea cluster, shows that the solution based on color features alone

AC

8

25

(a)

AN US

(b)

CR IP T

ACCEPTED MANUSCRIPT

Figure 12: Detailed view of the clusters from Figures 11(b) and (d) (k = 3). (a) Color features. (b) Color+texture features.

confuses part of the coastal line and terrain with the sea.

1

3.3. Profiling Meanshift

2

3

clusters as input. As a result, the STI profile is characterized by the number of

4

clusters as the independent variable, as shown in Sections 3.1 and 3.2. DBSCAN

5

ED

M

A number of clustering algorithms, such as k-means, require the number of

6

(Comaniciu and Meer , 2002) are examples of clustering algorithms that do not

7

require the number of clusters as input. Meanshift, in particular, was proposed

8

by Comaniciu and Meer (2002) as a clustering algorithm suitable for image

9

segmentation. Thus, we find it appropriate to analyze the STI profiles when

10

CE

PT

(Ester et al. , 1996), Affinity propagation (Frey and Dueck , 2007) and Meanshift

k-means is changed for Meanshift.

11

12

scale in the STI profile has real numbers pointing to the kernel bandwidths,

13

whereas the y-axis has the STI values corresponding to those bandwidths. The

14

new experiments are conducted on the natural images set. We aim at finding

15

the best bandwidth value and window size, hence, CS profiles need a normalized

16

x-axis calibration. Then, rather than using bandwidth values directly we use

17

percentiles from a distance distribution. Therefore, for each set of texture-

18

AC

Meanshift requires a kernel bandwidth as input parameter. Thus, the x-axis

26

CR IP T

ACCEPTED MANUSCRIPT

(b)

(c)

AN US

(a)

(d)

Figure 13: STI profiles for the color and color+texture experiments. (a) Peppers. (b) Limes.

M

(c) Landscape A. (d) Landscape B.

color and color features we calculate a distance matrix and construct a distance

2

histogram. As a result, instead of referring to a particular bandwidth value, we

3

can refer to the bandwidth value at percentile q. Histograms used as described

4

before normalize parameters within [0, 1], which makes it easier to compare

5

values that might have different scales.

PT

ED

1

6

Figure 14 shows the normalized bandwidth kernel STI profiles for two ex-

7

amples, Limes and Peppers. Normalized kernel bandwidth drives the number of clusters indirectly, i.e., the relation between number of clusters and bandwidth

9

values is not known a priori. However, we can rely on the STI index to pick

10

the optimum bandwidth and then inspect the results visually. Figures 15 and

AC

CE 8

11

16 are the gray scale masks that represent the best clustering solutions. These

12

masks were selected using the profiles displayed in Figure 14(a) and (c), which

13

describe the STI evolution for a bandwidth range. To improve visualization we

14

have included Figures 14(b) and (d), which magnify the portion of the profile

15

used to select the best Meanshift partitions. 27

CR IP T

ACCEPTED MANUSCRIPT

(b)

AN US

(a)

(c)

(d)

Figure 14: STI profiles for the color and color+texture experiments using Meanshift. (a)

M

Peppers. (b) Zoom in (a). (c) Limes. (d) Zoom in (c).

1

color and color+texture features. We highlighted the mask anomalies with red

2

ED

Segmentation masks from Figure 15 show similar partitions under different

3

number of anomalies, Figure 15(c), presents a slightly elevated STI index in the

4

profile displayed in Figure 14(b). In terms of STI values, Meanshift and k-means

5

PT

arrows. Interestingly, the partition corresponding to the mask with the most

6

the corresponding profiles are not correlated. Finally, a comparison between

7

the masks from Figure 7 (k-means) and Figure 15 (Meanshift) shows the two

8

different clustering algorithms arriving to almost identical results.

9

CE

best partitions are very similar. However, as Figures 13(a) and 14(a) show,

10

using only color features, Figure 16(a), has 4 clusters but the masks using

11

color+texture features, Figures 16(b-e), have 5 clusters. As in the previous

12

example, we have highlighted the most remarkable areas using red arrows. In

13

Figures 16(b) and (e) examples, the most salient distinction is an extra cluster

14

coming from the division of the Limes shadow. We point out this difference

15

AC

We display the masks dividing the Limes image in Figure 16. The mask

28

(b)

M

(a)

AN US

CR IP T

ACCEPTED MANUSCRIPT

(d)

PT

ED

(c)

(e)

CE

Figure 15: Peppers masks obtained using Meanshift. (a) Color features. (b) w = 31. (c)

AC

w = 41. (d) w = 51. (e) w = 61.

29

(b)

M

(a)

AN US

CR IP T

ACCEPTED MANUSCRIPT

(d)

PT

ED

(c)

(e)

CE

Figure 16: Limes masks obtained using Meanshift. (a) Color features. (b) w = 31. (c) w = 41.

AC

(d) w = 51. (e) w = 61.

30

ACCEPTED MANUSCRIPT

using (1) in Figure 16(b) and (2) in Figure 16(e). Figures 16(c) and (d), on

2

the other hand, have a small extra white cluster shown in (1) in both figures.

3

The best solutions achieved by CS and k-means are displayed in Figures 8 and

4

13(b), which select partitions with 4 clusters. These solutions are similar to

5

the one shown in Figure 16(a). However, the best Meanshift partition, which

6

appears at (0.2, w = 61) has a STI index similar to partitions using k-means

7

and 4 clusters.

8

3.4. Comparison with other state-of-the-art methods on a benchmark dataset

CR IP T

1

In this Section we compare CS combined with K-means as clustering algo-

10

rithm against GSWB (Mur et al. , 2016), Meanshift (Comaniciu and Meer ,

11

2002) and a new Level Sets based method, named Locally Statistical Active

12

Contour Model (LSACM) (Zhang et al. , 2016). For this purpose we use 4

13

gray scale natural images and 4 color natural images from the Berkeley Image

14

Dataset (Martin , 2001). The reason to select these images is they are the same

15

used by Mur et al.

16

We compare the performance of the four methods by visual inspection. Also,

17

for CS, GSWB and Meanshift, we measure the Silhouette index (Kaufman and

20

M

ED

19

(2016), whose method we use for comparison purposes.

Rousseeuw , 1990) and the accuracy value. The Silhouette index can be defined P as: Sk = N1 i s(i), where s(i) is the silhouette of point i, N is the number of patterns in the dataset and k the number of clusters in the solution. Next, we

PT

18

AN US

9

b(i)−a(i) max(b(i),a(i)) ,

21

define a single point i silhouette as: s(i) =

where a(i) represents

22

the average dissimilarity of i to all other patterns within the same cluster and b(i) is the smallest average dissimilarity of i to all the patterns of any other

24

cluster, of which i does not belong. Thus, s(i) belongs to: −1 < s(i) < 1, where

25

values near 1 represent good clustering and -1, bad clustering solutions. Like-

AC

CE 23

26

wise, the Silhouette index ranges from -1 to 1. The value interpretation is global

27

in this case rather than for a single point. Therefore, a Silhouette index value

28

of 1 indicates good clustering and a values of -1 indicate poor clustering. The

29

Silhouette index measures the clustering quality in the RGB or gray scale space

30

of the image, while the accuracy measures the percentage of pixels correctly 31

ACCEPTED MANUSCRIPT

classified. We used the manually segmentation provided by the Berkeley Image

1

Dataset to calculate the accuracy. When more than one solution was available,

2

we calculated the mean accuracy.

3

4

CR IP T

GSWB is a new method from a family of algorithms that uses spectral

5

This method automatically outputs a clustering solution that maximizes the

6

Silhouette index but using the transformed spectral space. We followed the

7

configuration proposed by Mur et al.

(2016) for GSWB with np = 700 pix-

8

els. On the other hand, Meanshift outputs a clustering solution for a given

9

bandwidth value. As to find a meaningful clustering solution, we used a grid

10

search over the bandwidth space [0.1, 1] while measuring the Silhouette index

11

to find the best clustering / segmentation. In CS case, we replicated the ex-

12

periment from Section 3.2 using the following window parameters for color and

13

color+texture features: w = {5, 11, 21, 31, 41, 51, color}. From the multiple con-

14

figurations, we kept a single solution with the minimum STI while measuring its

15

Silhouette index as well. Following the advice of Mur et al. (2016), we consid-

16

ered a simplified version of the Silhouette index as to speed up computations.

17

The accuracy values were measured over the same solutions as the Silhouettes,

18

ED

M

AN US

decomposition and intelligent image sampling to perform segmentation tasks.

19

we also compared our results to LSACM (Zhang et al. , 2016). Zhang et al.

20

(2016) successfully used their method to segment medical and natural images.

21

We show these results separately from the clustering methods since we only

22

compare its results visually.

23

CE

PT

which were contrasted against the ground truth provided by humans. Finally,

24

gray scale image experiments, GSWB and CS performed similarly, for instance

25

Image 1 and 2. For the rest of the images, many differences arise. First, in Image

26

3 CS divides the background into grass and flowers while GSWB segments the

27

background as a whole. Secongd, in Image 4, CS joins the church dome, the bells

28

and the windows with the background. On the other hand, GSWB divides part

29

of the dome and bells area into a separate cluster from the background. Finally,

30

in the Meanshift case, we can see that with the exception of Image 3 the rest

31

AC

After a visual inspection of the results, we can perceive that in some of the

32

(b) CS K-means

(e) Image 2

(f) CS K-means

(c) GSWB

(h) Meanshift

(j) CS K-means

(k) GSWB

(l) Meanshift

(n) CS K-means

(o) GSWB

(p) Meanshift

ED

M

(g) GSWB

CE

PT

(i) Image 3

(d) Meanshift

AN US

(a) Image 1

CR IP T

ACCEPTED MANUSCRIPT

(m) Image 4

AC

Figure 17: Segmentation results on four Berkeley gray scale images. Each row presents original image, CS K-means, GSWB and Meanshift results.

33

ACCEPTED MANUSCRIPT

Table 2: Silhouette indexes from CS K-means, GSWB and Meanshift results corresponding to the gray scale benchmark images from Figure 17. The number of clusters found by each method is reported between parentheses. We highlight the best result with Boldface fonts.

CS K-means

GSWB

Meanshift

Image 1

0.865 (3)

0.876 (3)

0.901 (2)

Image 2

0.694 (3)

0.729 (3)

0.707 (2)

Image 3

0.621 (3)

0.621 (2)

0.625 (2)

Image 4

0.864 (2)

0.784 (4)

0.851 (2)

AN US

CR IP T

Image

Table 3: Accuracy values from CS K-means, GSWB and Meanshift results corresponding to the gray scale benchmark images from Figure 17. Accuracies are in the range [0, 1]. We highlight the best result with Boldface fonts.

CS K-means

GSWB

Meanshift

Image 1

0.755

0.778

0.737

Image 2

0.661

0.659

0.567

0.649

0.718

0.661

0.637

0.652

0.631

Image 3

ED

Image 4

M

Image

PT

of the segmentations seem to be quite poor. Table 2 summarizes the Silhouette

1

2

in each test. Opposite to the perceptual analysis the Silhouette index measures

3

the clustering quality in terms of within cluster tightness and between cluster

4

separation. Table 2 shows that half of the examples are “better” clustered, in

5

terms of the Silhouette, by Meanshift. This could be expected since we used the

6

AC

CE

indexes for each method and image. It also reports the number of clusters found

Silhouette to find the best partition. The best solution for the two remaining

7

images are found by GSWB and CS. The accuracy values from Table 3, on the

8

contrary, show that GSWB performs better than CS in gray scale images.

9

We ran a second batch of simulations using a color version of the images

10

shown in Figure 17. From a perceptual point of view, CS shows a better per-

11

34

ACCEPTED MANUSCRIPT

Table 4: Silhouette indexes from CS K-means, GSWB and Meanshift results corresponding to the color benchmark images from Figure 18. The number of clusters found by each method is reported between parentheses. We highlight the best result with Boldface fonts.

CS K-means

GSWB

Meanshift

Image 5

0.855 (3)

0.828 (3)

0.891 (2)

Image 6

0.635 (3)

0.678 (3)

0.688 (2)

Image 7

0.312 (3)

0.535 (2)

0.483 (2)

Image 8

0.771 (4)

0.793 (3)

0.793 (4)

AN US

CR IP T

Image

Table 5: Accuracy values from CS K-means, GSWB and Meanshift results corresponding to the color benchmark images from Figure 18. Accuracies are in the range [0, 1]. We highlight the best result with Boldface fonts.

CS K-means

GSWB

Meanshift

Image 5

0.821

0.784

0.775

Image 6

0.744

0.715

0.646

0.720

0.667

0.637

0.748

0.740

0.747

Image 7

ED

Image 8

M

Image

formance than the other two methods. CS segmentation shows finer details and

2

sometimes finds objects not detected by either GSWB or Meanshift. For exam-

3

ple, CS Image 5 solution includes the moon while the other methods do not.

4

Also, CS Image 6 result includes the complete railing in the background and

5

the face details are more accurate. Moreover, CS Image 7 segmentation finds

6

both horses and divides grass and flowers. Finally, CS and Meanshift Image 8

CE

PT

1

results are quite similar but differ to GSWB result, which recognizes 3 clusters

8

rather than 4. Table 4 shows the Silhouette indexes and number of clusters

9

found for each clustering algorithm and image. Again, Meanshift outperforms

10

GSWB and CS in terms of Silhouette index. However, in this set of simulations,

11

CS presents the highest level of accuracy with respect to the ground truth.

AC 7

35

(b) CS K-means

(e) Image 6

(f) CS K-means

(c) GSWB

(d) Meanshift

AN US

(a) Image 5

CR IP T

ACCEPTED MANUSCRIPT

(h) Meanshift

(i) Image 7

(j) CS K-means

(k) GSWB

(l) Meanshift

PT

ED

M

(g) GSWB

(n) CS K-means

(o) GSWB

(p) Meanshift

CE

(m) Image 8

Figure 18: Segmentation results on four Berkeley color images. These images are the color

AC

version from the ones in Figure 17. Each row presents original image, CS K-means, GSWB and Meanshift results.

36

ACCEPTED MANUSCRIPT

(b)

(c)

AN US

CR IP T

(a)

(d)

Figure 19: LSACM (Zhang et al. , 2016) segmentation results for the four Berkeley gray scale images shown in Figure 17.

We also report the results obtained by the recently proposed Locally Statis-

2

tical Active Contour Model (LSACM) (Zhang et al. , 2016), which reported to

3

succesfully segment both medical and natural images. However, we only hoped

4

for partial comparison since this method can only process gray scale images.

5

Figure 19 shows the results on Images 1 to 4 after 1000 iterations. The method

6

was unable to successfully divide foreground from background in some of the

7

considered images. In the Image 2 example an its segmentation result Figure

8

19(c) we can find some similarities with the CS result depicted in Figure 17(h).

CE

PT

ED

M

1

Also, we find resemblance between Figure 19(d) (the segmentation result from

10

Image 3) and Figures 17(o) and 17(p), the results from GSBW and Meanshift,

11

respectively.

AC

9

12

3.5. Synthetic MRI

13

For the last set of experiments, we use k-means and Meanshift to parti-

14

tion/segment 3D synthetic MRI brain images. Our goal is to compare the per-

15

formance of clustering validation methods, nevertheless, the election of these 37

ACCEPTED MANUSCRIPT

1

(Tibshirani et al. , 2003) and Silhouette (Kaufman and Rousseeuw , 1990) are

2

computationally expensive both in memory and processing time. Thus, we se-

3

lected less complex methods such as WB (Zhao et al. , 2009) and CH (Calinski

4

CR IP T

methods is not unrestricted. Commonly used algorithms like Gap Statistic

and Harabasz , 1974) indexes to test against CS. Computational cost is an issue

5

for these kind of images since our examples have over 2 million voxels.

6

7

CH to find the number of clusters. Afterwards, we choose a winning partition

8

selected by each method to label the image voxels. To this end, we used the

9

same configuration of CS than before. Yet, for the WB and CH indexes we

10

use the procedure described by Zhao et al. (2009). To define WB and CH we

11

AN US

The experiments consist in clustering the MRI data using CS, WB and

first need to introduce the sum-of-squares within cluster (SSW) and the sum-ofP squares between clusters (SSB) metrics: SSWk = k,xi ∈Ck kxi − xk k2 , where xk is the center of cluster Ck , xi ∈ Ck and k = {1, 2, ..., K} is the number of P clusters for which we are calculating SSW. SSBk = k nk kxk − Xk2 , where xk

12

13

14

15

16

and k meaning is analogous to the one used in the previous formula. We now

17

define CH as: CHk =

M

is cluster k center, nk is cluster k cardinality, X represents the dataset mean SSBk /k−1 SSWk /N −k ,

18

we are calculating CH and N is the number of patterns in the dataset. Finally,

19

k we define the WB index as: W Bk = k SSW SSBk , where k is the number of clusters

20

for which we are calculating WB. In short, the minimum WB index points to

21

the best clustering partition. Hence, we can find this partition by measuring

22

the WB index within a range of (k = {2, 3, ..., 10}). On the contrary, the

23

maximum CH index indicates the best partitioning within a range of k. In the

24

Meanshift case, we measure the indexes within a bandwidth range instead, and

25

then we select the partition corresponding to the best index. Finally, we use

26

these three methods to select the best partitions according to each method, and

27

later, we compare these partitions against the ground truth. Our results report

28

the clustering accuracy and the number of clusters detected.

29

AC

CE

PT

ED

where k is the number of cluster for which

Table 6 shows the accuracy and number of clusters found by k-means when

30

using CS, WB and CH to find the best solution. Similarly, Table 7 presents

31

38

ACCEPTED MANUSCRIPT

Table 6: Accuracy and number of clusters for 3D MRI images. K-means+{CS, WB, CH}. The estimated number of clusters is shown in brackets.

CS

WB

BR 01

5

0.64 (6)

0.42 (10)

BR 02

5

0.83 (4)

0.41 (10)

BR 03

5

0.73 (5)

0.41 (10)

BR 04

5

0.70 (5)

0.40 (10)

BR 05

5

0.72 (5)

0.42 (10)

BR 06

5

0.67 (5)

0.37 (10)

CH

CR IP T

# Classes

0.40 (10) 0.42 (10) 0.41 (10) 0.40 (10) 0.42 (10) 0.37 (10)

AN US

Image

the same information but using Meanshift as clustering algorithm. Unlike our

2

previous test with real images, simulated MRI data is generated with ground

3

truth. Thus, we can evaluate the performance of any method directly, without

4

relying on visual perception. The k-means partitions selected by WB and CH

5

indexes (see Table 6) tend to fragment the solution into many spurious clusters.

6

On the other hand, the number of clusters found by CS are more representative

7

of the actual number of classes. From the accuracy levels we infer that the

8

partitions selected by the WB and CH indexes are the outcome of the true

9

classes being heavily fragmented. On the contrary, if extra clusters containing

10

small numbers of voxels were detected, the accuracy levels would be higher.

11

Accuracy levels from Table 6 indicate that partitions selected by CS are more

12

correlated to true classes than solutions found by WB an CH.

CE

PT

ED

M

1

Table 7 shows the accuracy based on Meanshift partitions validated by CS,

14

WB and CH. Most of the cases indicate that Meanshift solutions tend to un-

15

derestimate the number of classes. Anyways, accuracies are better than using

16

k-means. The reason for having higher accuracies, in spite of underestimat-

17

ing the number of classes, is because classes are highly imbalanced, voxels are

18

mostly distributed between: white matter, gray matter and cerebrospinal fluid.

19

Thus, detecting these three classes already ensures high accuracy.

AC

13

39

ACCEPTED MANUSCRIPT

Table 7: Accuracy and number of clusters for 3D MRI images. Meanshift+{CS, WB, CH}. The estimated number of clusters is shown in brackets.

CS

WB

BR 01

5

0.79 (4)

0.91 (4)

BR 02

5

0.89 (3)

0.62 (8)

BR 03

5

0.90 (3)

0.83 (5)

BR 04

5

0.81 (4)

0.88 (4)

BR 05

5

0.88 (3)

0.83 (5)

BR 06

5

0.81 (3)

0.37 (9)

CH

CR IP T

# Classes

0.91 (4) 0.62 (8) 0.90 (4) 0.88 (4) 0.89 (4) 0.81 (4)

AN US

Image

1

database, as well as the segmentation results obtained by CS, WB and CH using

2

k-means and Meanshift (Figures 20(b-g)). These results confirm the previous

3

discussion about the performance achieved by of each one of the compared

4

methods.

5

M

Figure 20(a) shows a typical brain image slice from the synthetic 3D MRI

6

pairwise distances to find the within cluster sum of squares (WSS) and between

7

clusters sum of squares (BSS), which requires quadratic time respect to the

8

number of points per cluster. Meanwhile, CS complexity just increases linearly

9

with the number of repetitions. But it also depends on the complexity of the

10

clustering algorithm and classifier selected. In these set of examples CS required

11

approximately 5 to 7 minutes while WB and CH needed at least 3 hours.

12

CE

PT

ED

An important disadvantage of WB and CH indexes is the need to calculate

AC

4. Concluding remarks and future work

13

Clustering is a well-established technique for segmentation. However, clus-

14

tering validation is rarely used either to support a solution based on clustering

15

or to automate the segmentation process. Therefore, our work is a step forward

16

to improve segmentation by clustering. The proposed method, not only intro-

17

duces some major modifications to CS, but also discusses a new representation

18

40

M

AN US

(a)

CR IP T

ACCEPTED MANUSCRIPT

(c)

(d)

(f)

(g)

CE

PT

ED

(b)

AC

(e)

Figure 20: Segmentation results for a typical synthetic MRI image using CS, WB and CH in combination with k-means and Meanshift. (a) Brain axial slice. (b) CS+k-means. (c) WB+kmeans. (d) CH+k-means. (e) CS+Meanshift. (f) WB+Meanshift. (g) CH+Meanshift.

41

ACCEPTED MANUSCRIPT

1

on texture and natural images, we were able to automatically find the “best”

2

image segmentation and texture parameters simultaneously, where “best” refers

3

to a solution that can be repeatedly learned by a classifier with the lowest er-

4

CR IP T

of the color+texture features to achieve memory efficiency. In our experiments

ror and it is not used as an absolut term. In the additional file, we included

5

a Gabor texture descriptor to expand our experimental results. These results

6

show that CS is not limited to the NRI-LBP texture features. Moreover, the

7

CS framework can select meaningful clustering results using the Gabor texture

8

descriptor as well. In like manner, other texture descriptors can be combined

9

AN US

with RGB information within the CS framework.

10

11

sults. For gray scale and color images Meanshift achieves the “best” clustering

12

according to the Silhouette index. In contrast, GSWB finds the clustering so-

13

lution with highest accuracy for gray scale images, followed by CS. Finally, on

14

color image experiments the highest accuracy was obtained by CS. The two best

15

methods, GSWB and CS, use a non-trivial validation procedure embedded with

16

the clustering rather than measuring a validation index directly to the data.

17

Results show that this strategy is more successful.

18

ED

M

The experiments performed on Berkeley Dataset images show interesting re-

19

be used to point out clutering solutions. We have presented experiments using

20

k-means, which was selected due to its speed. Interestingly, a second clustering

21

algorithm, Meanshift, achieved similar segmentation results. Finally, we ran

22

tests on 3D artificial MRI data. We used ground truth labels to compare the

23

CE

PT

CS uses a clustering algorithm and a classifier to find the STI, which can

24

CS provided a fast and accurate alternative. However, we note that clustering

25

algorithms had more influence on the results than in previous experiments.

26

AC

segmentation quality obtained by three clustering validation methods, where

A motivating challenge is to accurately achieve fully automated MRI data

27

segmentation by means of an unsupervised method such as CS. We are looking

28

into: i) improving data representation and ii) fusing several MRI modalities

29

(T1, T2, etc.) to improve segmentation.

30

One of the most interesting features of CS is that it is not limited by dimen42

31

ACCEPTED MANUSCRIPT

sionality, hence, 2D and 3D images can be exchanged easily. Additionally, any

2

clustering algorithm can be used. Moreover, if several algorithms are used, the

3

STI index allows to compare them and select the best in the same manner as

4

we did for the window sizes.

5

Acknowledgements

CR IP T

1

Authors acknowledge grant support from ANPCyT PICT-2012-0181.

6

References

8

Alata, O., Ramananjarasoa, C.: Unsupervised textured image segmentation

9

using 2-D quarter plane autoregressive model with four prediction supports.

AN US

7

Pattern Recognition Letters 26(8), 1069–1081 (2005).

10

Altuvia, Y., Landgraf, P., Lithwick, G., Elefant, N., Pfeffer, S., Aravin, A.,

12

Brownstein, M., Tuschl, T., Margalit, H.: Clustering and conservation pat-

13

terns of human microRNAs. Nucleic Acids Research 33(8), 2697–2706 (2005)

14

Arbel´ aez, P., Maire, M., Fowlkes, C., Malik, J.: Contour detection and hi-

15

erarchical image segmentation. IEEE Transactions on Pattern Analysis and

16

Machine Intelligence 33(5), 898–916 (2011).

ED

Aue, A., Lee, T.: On image segmentation using information theoretic criteria.

PT

17

M

11

Annals of Statistics 39(6), 2912–2935 (2011).

18

Bay´ a, A., Granitto, P.: How many clusters: A validation index for arbitrary-

CE

19

20

shaped clusters. IEEE/ACM Transactions on Computational Biology and

21

Bioinformatics 10(2), 401–414 (2013).

Bergasa, L., Mazo, M., Gardel, A., Sotelo, M., Boquete, L.: Unsupervised and

23

adaptive Gaussian skin-color model. Image and Vision Computing 18(12),

24

987–1003 (2000).

AC 22

25

26

Calinski, T., Harabasz, J.: A Dendrite Method for Cluster Analysis. Communications in Statistics 3, 1–27 (1974) 43

ACCEPTED MANUSCRIPT

Chen, J., Pappas, T., Mojsilovic, A., Rogowitz, B.: Adaptive perceptual color-

1

texture image segmentation. IEEE Transactions on Image Processing 14(10),

2

1524–1536 (2005)

3

4

Web: Online Interface to a 3D MRI Simulated Brain Database. NeuroImage

5

5, 425 (1997)

6

CR IP T

Cocosco, C.A., Kollokian, V., Kwan, R.K.S., Pike, G.B., Evans, A.C.: Brain-

7

analysis. IEEE Transactions on Pattern Analysis and Machine Intelligence

8

24(5), 603–619 (2002).

9

AN US

Comaniciu, D., Meer, P.: Mean shift: A robust approach toward feature space

de Luis-Garc´ıa, R., Deriche, R., Alberola-L´ opez, C.: Texture and color seg-

10

mentation based on the combined use of the structure tensor and the image

11

components. Signal Processing 88(4), 776–795 (2008).

12

Dhillon, I., Modha, D.: Concept decompositions for large sparse text data using

M

clustering. Machine Learning 42(1-2), 143–175 (2001)

13

14

15

discovering clusters in large spatial databases with noise. In: , pp. 226–231.

16

ED

Ester, M., Kriegel, H.P., Sander, J., Xu, X.: A density-based algorithm for

AAAI Press (1996)

17

Felzenszwalb, P., Huttenlocher, D.: Efficient graph-based image segmentation.

PT

International Journal of Computer Vision 59(2), 167–181 (2004). Frey, B.J., Dueck, D.: Clustering by Passing Messages Between Data Points.

CE

Science 315(5814), 972–976 (2007).

19

20

21

Handl, J., Knowles, J., Kell, D.B.: Computational cluster validation in post-

AC

18

genomic data analysis. Bioinformatics 21(15), 3201–3212 (2005)

Hastie, T., Tibshirani, R., Friedman, J.: The Elements of Statistical Learning: Data Mining, Inference, and Prediction. Springer (2003)

Herbin, M., Bonnet, N., Vautrot, P.: Estimation of the number of clusters and influence zones. Pattern Recognition Letters 22(14), 1557–1568 (2001). 44

22

23

24

25

26

27

ACCEPTED MANUSCRIPT

1

Ilea, D.E., Whelan, P.F.: CTex - An adaptive unsupervised segmentation al-

2

gorithm based on color-texture coherence. IEEE Transactions on Image Pro-

3

cessing 17(10), 1926–1939 (2008). Ilea, D.E., Whelan, P.F.: Image segmentation based on the integration of colour-

5

texture descriptors - A review. Pattern Recognition 44(10-11), 2479–2501

6

(2011).

CR IP T

4

Jakobsson, M., Rosenberg, N.: CLUMPP: A cluster matching and permutation

8

program for dealing with label switching and multimodality in analysis of

9

population structure. Bioinformatics 23(14), 1801–1806 (2007)

10

AN US

7

Kaufman, L., Rousseeuw, P.J.: Finding Groups in Data: An Introduction to Cluster Analysis. Wiley-Interscience (1990)

11

Khanmohammadi, S., Adibeig, N., Shanehbandy, S.: An improved overlapping

13

k-means clustering method for medical applications. Expert Systems with

14

Applications 67, 12–18 (2017)

15

M

12

Kim, J.-S., Hong, K.-S.: Color-texture segmentation using unsupervised graph cuts. Pattern Recognition 42(5), 735–750 (2009).

17

ED

16

Kuhn, H.W.: The Hungarian Method for the Assignment Problem. Naval Research Logistics Quarterly pp. 83–97 (1955)

19

PT

18

Lange, T., Roth, V., Braun, M.L., Buhmann, J.M.: Stability-Based Validation of Clustering Solutions. Neural Comp. 16(6), 1299–1323 (2004)

CE

20

Macqueen, J.: Some methods for classification and analysis of multivariate ob-

22

servations. In: In 5-th Berkeley Symposium on Mathematical Statistics and

AC

21

23

Probability, pp. 281–297 (1967)

24

Martin D., Fowlkes C., Tal D., Malik J.: A Database of Human Segmented Nat-

25

ural Images and its Application to Evaluating Segmentation Algorithms and

26

Measuring Ecological Statistics. In: Proc. 8th Int0l Conf. Computer Vision

27

Vol.(2), pp. 416–423 (2001) 45

ACCEPTED MANUSCRIPT

Mur, A., Dormido, R., Duro, N., Dormido-Canto, S., Vega, J.: Determination

1

of the optimal number of clusters using a spectral clustering optimization.

2

Expert Systems with Applications 65, 304–314 (2016)

3

4

relationship management: A literature review and classification. Expert Sys-

5

tems with Applications 36(2), 2592–2602 (2009).

6

CR IP T

Ngai, E., Xiu, L., Chau, D.: Application of data mining techniques in customer

7

tion invariant texture classification with local binary patterns. IEEE Trans-

8

actions on Pattern Analysis and Machine Intelligence 24(7), 971–987 (2002).

9

AN US

Ojala, T., Pietik¨ ainen, M., M¨ aenp¨ a¨ a, T.: Multiresolution gray-scale and rota-

Rousseeuw, P.: Silhouettes: a graphical aid to the interpretation and validation of cluster analysis. J. Comput. Appl. Math. 20, 53–65 (1987)

Shi, J., Malik, J.: Normalized cuts and image segmentation. IEEE Transactions on Pattern Analysis and Machine Intelligence 22(8), 888–905 (2000).

10

11

12

13

14

covery and Applications of Usage Patterns from Web Data. SIGKDD Explor.

15

Newsl. 1(2), 12–23 (2000).

16

ED

M

Srivastava, J., Cooley, R., Deshpande, M., Tan, P.N.: Web Usage Mining: Dis-

17

The pseudolikelihood information criterion (PLIC). IEEE Transactions on

18

PT

Stanford, D., Raftery, A.: Approximate bayes factors for image segmentation:

Pattern Analysis and Machine Intelligence 24(11), 1517–1520 (2002)

CE

Tarabalka, Y., Benediktsson, J., Chanussot, J.: Spectral-spatial classification

19

20

21

Transactions on Geoscience and Remote Sensing 47(8), 2973–2987 (2009).

22

Tibshirani, R., Walther, G., Hastie, T.: Estimating the Number of Clusters in

23

a Dataset Via the Gap Statistic. Journal of the Royal Statistical Society B

24

63, 411–423 (2003)

25

AC

of hyperspectral imagery based on partitional clustering techniques. IEEE

Tseng, Y.H., Lin, C.J., Lin, Y.I.: Text mining techniques for patent analysis. Information Processing and Management 43(5), 1216–1247 (2007) 46

26

27

ACCEPTED MANUSCRIPT

Vinh, N.X., Epps, J., Bailey, J.: Information Theoretic Measures for Clusterings

2

Comparison: Is a Correction for Chance Necessary? In: Proceedings of the

3

26th Annual International Conference on Machine Learning, ICML ’09, pp.

4

1073–1080 (2009).

CR IP T

1

5

Yang, Y., Han, S., Wang, T., Tao, W., Tai, X.C.: Multilayer graph cuts based

6

unsupervised color–texture image segmentation using multivariate mixed stu-

7

dent’s t-distribution and regional credibility merging. Pattern Recognition

8

46(4), 1101 – 1124 (2013).

Yeung, K.Y., Haynor, D.R., Ruzzo, W.L.: Validating Clustering for Gene Ex-

AN US

9

pression Data. Bioinformatics 17(Suppl 4) (2001)

10

11

Zhang K., Zhang L., Lam K. M., Zhang D.: A Level Set Approach to Image

12

Segmentation With Intensity Inhomogeneity. IEEE Transactions on Cyber-

13

netics,46(2), 546–557 (2016).

Zhao, Q., Xu, M., Fr¨ anti, P.: Sum-of-squares based cluster validity index and

15

significance analysis. Lecture Notes in Computer Science (including subseries

16

Lecture Notes in Artificial Intelligence and Lecture Notes in Bioinformatics)

17

5495 LNCS, 313–322 (2009).

AC

CE

PT

ED

M

14

47