Clustering an image using Gaussian mixture models

Question

I want to use GMM(Gaussian mixture models for clustering a binary image and also want to plot the cluster centroids on the binary image itself.

I am using this as my reference: http://in.mathworks.com/help/stats/gaussian-mixture-models.html

This is my initial code

 I=im2double(imread('sil10001.pbm'));
  K = I(:);
  mu=mean(K);
  sigma=std(K);
  P=normpdf(K, mu, sigma);
   Z = norminv(P,mu,sigma);
  X = mvnrnd(mu,sigma,1110);
  X=reshape(X,111,10);


 scatter(X(:,1),X(:,2),10,'ko');

options = statset('Display','final');
gm = fitgmdist(X,2,'Options',options);



idx = cluster(gm,X);
cluster1 = (idx == 1);
cluster2 = (idx == 2);


 scatter(X(cluster1,1),X(cluster1,2),10,'r+');
 hold on

  scatter(X(cluster2,1),X(cluster2,2),10,'bo');
  hold off
  legend('Cluster 1','Cluster 2','Location','NW')


  P = posterior(gm,X);

 scatter(X(cluster1,1),X(cluster1,2),10,P(cluster1,1),'+')
 hold on
 scatter(X(cluster2,1),X(cluster2,2),10,P(cluster2,1),'o')
 hold off
 legend('Cluster 1','Cluster 2','Location','NW')
 clrmap = jet(80); colormap(clrmap(9:72,:))
 ylabel(colorbar,'Component 1 Posterior Probability')

But the problem is that I am unable to plot the cluster centroids received from GMM in the primary binary image.How do i do this?

**Now suppose i have 10 such images in a sequence And i want to store the information of their mean position in two cell array then how do i do that.This is my code foe my new question **

    images=load('gait2go.mat');%load the matrix file
    for i=1:10

   I{i}=images.result{i};
  I{i}=im2double(I{i});

   %determine 'white' pixels, size of image can be [M N], [M N 3] or [M N 4]
  Idims=size(I{i});
  whites=true(Idims(1),Idims(2));

    df=I{i};
      %we add up the various color channels
 for colori=1:size(df,3)
  whites=whites & df(:,:,colori)>0.5;
 end

%choose indices of 'white' pixels as coordinates of data
[datax datay]=find(whites);

%cluster data into 10 clumps
  K = 10;               % number of mixtures/clusters
  cInd = kmeans([datax datay], K, 'EmptyAction','singleton',...
   'maxiter',1000,'start','cluster');

%get clusterwise means
 meanx=zeros(K,1);
 meany=zeros(K,1);  
  for i=1:K
   meanx(i)=mean(datax(cInd==i));
   meany(i)=mean(datay(cInd==i));

 end

 xc{i}=meanx(i);%cell array contaning the position of the mean for the 10    
 images
  xb{i}=meany(i);

figure;
gscatter(datay,-datax,cInd); %funky coordinates for plotting according to      
 image
 axis equal;
  hold on;
  scatter(meany,-meanx,20,'+'); %same funky coordinates


 end

I am able to get 10 images segmented but no the values of themean stored in the cell arrays xc and xb.They r only storing [] in place of the values of means

Answer 1

I decided to post an answer to your question (where your question was determined by a maximum-likelihood guess:P), but I wrote an extensive introduction. Please read carefully, as I think you have difficulties understanding the methods you want to use, and you have difficulties understanding why others can't help you with your usual approach of asking questions. There are several problems with your question, both code-related and conceptual. Let's start with the latter.

The problem with the problem

You say that you want to cluster your image with Gaussian mixture modelling. While I'm generally not familiar with clustering, after a look through your reference and the wonderful SO answer you cited elsewhere (and a quick 101 from @rayryeng) I think you are on the wrong track altogether.

Gaussian mixture modelling, as its name suggests, models your data set with a mixture of Gaussian (i.e. normal) distributions. The reason for the popularity of this method is that when you do measurements of all sorts of quantities, in many cases you will find that your data is mostly distributed like a normal distribution (which is actually the reason why it's called normal). The reason behind this is the central limit theorem, which implies that the sum of reasonably independent random variables tends to be normal in many cases.

Now, clustering, on the other hand, simply means separating your data set into disjoint smaller bunches based on some criteria. The main criterion is usually (some kind of) distance, so you want to find "close lumps of data" in your larger data set. You usually need to cluster your data before performing a GMM, because it's already hard enough to find the Gaussians underlying your data without having to guess the clusters too. I'm not familiar enough with the procedures involved to tell how well GMM algorithms can work if you just let them work on your raw data (but I expect that many implementations start with a clustering step anyway).

To get closer to your question: I guess you want to do some kind of image recognition. Looking at the picture, you want to get more strongly correlated lumps. This is clustering. If you look at a picture of a zoo, you'll see, say, an elephant and a snake. Both have their distinct shapes, and they are well separated from one another. If you cluster your image (and the snake is not riding the elephant, neither did it eat it), you'll find two lumps: one lump elephant-shaped, and one lump snake-shaped. Now, it wouldn't make sense to use GMM on these data sets: elephants, and especially snakes, are not shaped like multivariate Gaussian distributions. But you don't need this in the first place, if you just want to know where the distinct animals are located in your picture.

Still staying with the example, you should make sure that you cluster your data into an appropriate number of subsets. If you try to cluster your zoo picture into 3 clusters, you might get a second, spurious snake: the nose of the elephant. With an increasing number of clusters your partitioning might make less and less sense.

Your approach

Your code doesn't give you anything reasonable, and there's a very good reason for that: it doesn't make sense from the start. Look at the beginning:

I=im2double(imread('sil10001.pbm'));
K = I(:);
mu=mean(K);
sigma=std(K);
X = mvnrnd(mu,sigma,1110);
X=reshape(X,111,10);

You read your binary image, convert it to double, then stretch it out into a vector and compute the mean and deviation of that vector. You basically smear your intire image into 2 values: an average intensity and a deviation. And THEN you generate 111*10 standard normal points with these parameters, and try to do GMM on the first two sets of 111. Which are both independently normal with the same parameter. So you probably get two overlapping Gaussians around the same mean with the same deviation.

I think the examples you found online confused you. When you do GMM, you already have your data, so no pseudo-normal numbers should be involved. But when people post examples, they also try to provide reproducible inputs (well, some of them do, nudge nudge wink wink). A simple method for this is to generate a union of simple Gaussians, which can then be fed into GMM.

So, my point is, that you don't have to generate random numbers, but have to use the image data itself as input to your procedure. And you probably just want to cluster your image, instead of actually using GMM to draw potatoes over your cluster, since you want to cluster body parts in an image about a human. Most body parts are not shaped like multivariate Gaussians (with a few distinct exceptions for men and women).

What I think you should do

If you really want to cluster your image, like in the figure you added to your question, then you should use a method like k-means. But then again, you already have a program that does that, don't you? So I don't really think I can answer the question saying "How can I cluster my image with GMM?". Instead, here's an answer to "How can I cluster my image?" with k-means, but at least there will be a piece of code here.

%set infile to what your image file will be
infile='sil10001.pbm';

%read file
I=im2double(imread(infile));

%determine 'white' pixels, size of image can be [M N], [M N 3] or [M N 4]
Idims=size(I);
whites=true(Idims(1),Idims(2));

%we add up the various color channels
for colori=1:Idims(3)
    whites=whites & I(:,:,colori)>0.5;
end

%choose indices of 'white' pixels as coordinates of data
[datax datay]=find(whites);

%cluster data into 10 clumps
K = 10;               % number of mixtures/clusters
cInd = kmeans([datax datay], K, 'EmptyAction','singleton',...
    'maxiter',1000,'start','cluster');

%get clusterwise means
meanx=zeros(K,1);
meany=zeros(K,1);
for i=1:K
    meanx(i)=mean(datax(cInd==i));
    meany(i)=mean(datay(cInd==i));
end

figure;
gscatter(datay,-datax,cInd); %funky coordinates for plotting according to image
axis equal;
hold on;
scatter(meany,-meanx,20,'ko'); %same funky coordinates

Here's what this does. It first reads your image as double like yours did. Then it tries to determine "white" pixels by checking that each color channel (of which can be either 1, 3 or 4) is brighter than 0.5. Then your input data points to the clustering will be the x and y "coordinates" (i.e. indices) of your white pixels.

Next it does the clustering via kmeans. This part of the code is loosely based on the already cited answer of Amro. I had to set a large maximal number of iterations, as the problem is ill-posed in the sense that there aren't 10 clear clusters in the picture. Then we compute the mean for each cluster, and plot the clusters with gscatter, and the means with scatter. Note that in order to have the picture facing in the right directions in a scatter plot you have to shift around the input coordinates. Alternatively you could define datax and datay correspondingly at the beginning.

And here's my output, run with the already processed figure you provided in your question: