Sort clusters using K-means by intensity

Hello everyone. I am using K-means to segment some grayscale images. Unfortunately, the values of the generated clusters are not repeatable, i.e. every time I run the code the clusters have a different value. For example, if I use k=2 sometimes the darker areas of the original image have a cluster value of 1 and sometimes 2 (before normalisation). How to sort/order the generated clusters to have a value corresponding to the actual grayscale intensities, i.e. darkest = 1, less dark = 2,... brightest = k ? Thanks. Here is the code:

      % Clustering.
      clustered = reshape(kmeans(inputimage(:), k), size(inputimage)); 

      % Normalise intensities from 0 to 1.
      clustered = clustered - min(clustered(:));
      clustered = clustered / max(clustered(:));

[Xen]

Anybody???

Answer 1

[Walter Roberson]

You are normalizing the indices, not by cluster intensities.

kidx = kmeans(inputimage(:), k);
clustermeans = accumarray(kidx, inputimage(kidx),[], @mean);
[sortedmeans, sortidx] = sort(clustermeans);
kidxmapped = sortidx(kidx);
clustered = reshape(kidxmapped, size(inputImage));

[Xen]

Hi Walter. I am not sure why, but it didn't work. I still get mixed clusters.

[Walter Roberson]

kidx = kmeans(inputimage(:), k);
clustermeans = accumarray(kidx, inputimage(kidx),[], @mean);
[sortedmeans, ~, sortidx] = unique(clustermeans);
kidxmapped = sortidx(kidx);
clustered = reshape(kidxmapped, size(inputImage));

[Xen]

Thanks for trying to help, but still, neither this worked. The code seems fine and I am busting my head to figure out why it doesn't work. I am showing below the output from some of the functions, for k=3.

clustermeans =

      0.0104
      0.0063
      0.0005

sortedmeans =

      0.0005
      0.0063
      0.0104

sortidx =

       3
       2
       1

This is exactly the same every time I run the code, and while it looks it does the correct thing, it doesn't. For example, when I estimate the mode(clustered(:)) at the end I sometimes get 1, sometimes 2, sometimes 3, and the output image looks different accordingly (with exactly the same input image). I'm freaking out! Can you try this on an image and tell me whether I am the only one who get's this?

[Walter Roberson]

kmeans uses random numbers to initialize the cluster centers. There is an option to kmeans to provide initial cluster centers. If you want to use that, pass kmeans the 'Start' option followed by,

k-by-p matrix of centroid starting locations. The rows of Start correspond to seeds. The software infers k from the first dimension of Start, so you can pass in [] for k.

mode() on the cluster indices is essentially going to tell you which cluster had the most elements. kmeans is deterministic once the initial "seeds" are known, but if you do not specify the seeds then it is not deterministic as to where it will place them.

You might wish to look at the Replicates option, which tells kmeans to run several times with different seeds and choose the best outcome.

kmeans does not find a global minimum, it finds a local minimum, so the starting point is important.

[Xen]

I used mode() just to see in a faster way whether the code gives me the same image, without having to look at it. I understand that kmeans initializes at random centers, but regardless of this I would expect

kidxmapped = sortidx(kidx);

to sort things out and give me the correct image. I can't understand why it fails.

Thanks anyway. I'll see if another option can solve this.

[Walter Roberson]

I think you might be a lot happier using

thresh = multithresh(inputimage, k-1);
[counts, kidx] = histc(inputImage(:), [-inf thresh inf]);
clustermeans = accumarray(kidx(:), inputimage(kidx),[], @mean);
clustered = reshape(kidx, size(inputimage));

The clusters should already be in sorted order by mean.

The method used for multithresh is deterministic.

[Xen]

Multithreshold works, yes. Regarding kmeans (I am very persistent, I know), I've found out that sort() gives me inappropriate indexes of the sorted cluster means. For example,

clustermeans =

      0.1401
      0.2703
      0.0141

sortidx =

       3
       1
       2

However, sortidx must give me the indices of the original array as if it was expanded from 1 to 3, and not the indices of how to arrange clustermeans, i.e.

2
3
1

Is there a function to actually do this?

[Xen]

Yeeess! I finally cracked this. I created a method to sort the clustermeans appropriately, and also had to use a different method to create clustermeans otherwise it kept failing. This is quite slow for large k values. If you have any suggestions for improvement I would be glad to know.

    % Cluster.
    clustered = reshape(kmeans(image(:), k), size(image));

    % Sort clusters.
    clustermeans = zeros(k, 1);
    for j = 1:k
        [row, col] = find(clustered == j);
        pixelsKvalued = zeros(length(row), 1);
        for l = 1:length(row)
            pixelsKvalued(l) = image(row(l), col(l));
        end
        clustermeans(j) = mean(pixelsKvalued);
    end 

    sortidx = zeros(k, 1);
    for j = 1:k
        sortidx(find(clustermeans == min(clustermeans))) = j;    
        clustermeans(clustermeans == min(clustermeans)) = NaN;
    end

    clustered = sortidx(clustered);

[Xen]

No need to estimate the mean of all pixels representing each cluster. Just a single pixel will do. Thanks Walter, this was inspired by your suggestions.

    % Cluster.
    clustered = reshape(kmeans(image(:), k), size(image));

    % Sort clusters. 
    clusterintensity = zeros(k, 1);
    for j = 1:k
        clusterintensity(j) = image(find(clustered == j, 1));
    end 

    clusteridx = zeros(k, 1);
    for j = 1:k
        clusteridx(find(clusterintensity == min(clusterintensity))) = j;    
        clusterintensity(clusterintensity == min(clusterintensity)) = NaN;
    end

    clustered = clusteridx(clustered);