How to fit Gaussian mixture models on image histogram?

Hello all,

I am writing a code to seperate clusters of an image. I did this part. However, I want to show the GMM on the image Histogram. Here the code that I used:

% Read Image


I = im2double(imread('test.tif'));

mask = I > 0 & I < 65536; %remove background

I1 = I(mask);
data = I1(:);

% Fit a gaussian mixture model
n_components = 4;
obj = fitgmdist(data,n_components);
idx = cluster(obj,data);
cluster1 = data(idx == 1,:);
cluster2 = data(idx == 2,:);
cluster3 = data(idx == 3,:);
cluster4 = data(idx == 4,:);

% Display Histogram

histogram(cluster1,0:.01:1); hold on;
histogram(cluster2,0:.01:1); hold on;
histogram(cluster3,0:.01:1); hold on;
histogram(cluster4,0:.01:1);

I would also like to show each cluster as a seperate Gaussian curve.

Thanks,

Answer 1

[Akshat]

Hi Mohamed,

I understand that you would like to plot a depiction of the GMM over the image histogram you are generating in the code you mentioned. I assume that the PDF (Probability Density Function) is a good depiction of the GMM.

As you did not provide any image for the question, I generated an image using the code attached as image_gen.m. Kindly refer it in case you want to know what the image looks like.

Below is the modified code:

% Read Image 

I = im2double(imread('test.tif')); 

 

mask = I > 0 & I < 65536; % Remove background 

I1 = I(mask); 

data = I1(:); 

 

% Fit a Gaussian mixture model 

n_components = 4; 

obj = fitgmdist(data, n_components); 

idx = cluster(obj, data); 

cluster1 = data(idx == 1, :); 

cluster2 = data(idx == 2, :); 

cluster3 = data(idx == 3, :); 

cluster4 = data(idx == 4, :); 

 

% Display Histogram 

 

histogram(cluster1, 0:.01:1); hold on; 

histogram(cluster2, 0:.01:1); hold on; 

histogram(cluster3, 0:.01:1); hold on; 

histogram(cluster4, 0:.01:1); hold off; 

 

% Plot GMM PDFs on the histogram 

x = linspace(0, 1, 101); % Values for x-axis 

x = x';  

pdf1 = pdf(obj, x); % GMM PDF for cluster 1 

pdf2 = pdf(obj, x); % GMM PDF for cluster 2 

pdf3 = pdf(obj, x); % GMM PDF for cluster 3 

pdf4 = pdf(obj, x); % GMM PDF for cluster 4 

 

binEdges = 0:.01:1; 

figure; 

histogram(cluster1, binEdges); hold on; 

plot(x, pdf1 * numel(cluster1) * diff(binEdges(1:2)), 'r', 'LineWidth', 2); 

histogram(cluster2, binEdges); hold on; 

plot(x, pdf2 * numel(cluster2) * diff(binEdges(1:2)), 'g', 'LineWidth', 2); 

histogram(cluster3, binEdges); hold on; 

plot(x, pdf3 * numel(cluster3) * diff(binEdges(1:2)), 'b', 'LineWidth', 2); 

histogram(cluster4, binEdges); hold on; 

plot(x, pdf4 * numel(cluster4) * diff(binEdges(1:2)), 'm', 'LineWidth', 2); 

hold off; 

% Set plot labels and title 

xlabel('Pixel Value'); 

ylabel('Frequency'); 

title('Image Histogram with GMM PDFs'); 

legend('Cluster 1', 'Cluster 1 PDF', 'Cluster 2', 'Cluster 2 PDF', 'Cluster 3', 'Cluster 3 PDF', 'Cluster 4', 'Cluster 4 PDF'); 

 

Hope this helps!