Image segmentation by k-means algorithm

I have an code for k-means segmentation but I have some problems when I applying it on my multispectral satellite image I facing error in dimensions (Subscripted assignment dimension mismatch), how can I solve this problem?

%  Grayscale Image Segmentation Using K-Means Algorithm Function Kmeans segmentation

clc
close all
clear all

[im,map]=imread('pp1.bmp');
im=im2double(im);
[row,col]=size(im);

%  number of clusters
nc=4;

%  initial random cluster centroids
cs=rand(nc,1);
pcs=cs;

%  number of iteration 
T=50;
t=0;
D=zeros(row,col,nc);
tsmld=[];
eps=1.e-5;
cmx=1;

while (t<T && cmx>eps)
    %Distance between centroids and image's pixel
    for c=1:nc
        D(:,:,c)=(im-cs(c)).^2;
    end

      %assign members (image pixels)to minimum distance clusters
      [mv,ML]=min(D,[],3);

      %updat cluster centroid
      for c=1:nc
          I=(ML==c);
          cs(c)=mean(mean(im(I)));
      end

      %find maximum absolute difference between crrent and previous iteration
      %cluster centroids
      cmx=max(abs(cs-pcs));
      pcs=cs;

      t=t+1;

      %sum difference between centroid and their members and store it for
      %plotting energy minimization functions
      tsmld=[tsmld; sum(mv(:))];

end

%  assign a colour to each cluster

colors=hvs(nc);
sim=colors(ML,:);
sim=reshape(sim,row,col,3);

figure,subplot(1,2,1),imshow(im,map);
title('Input Image: pp1');
subplot(1,2,2);imshow(sim,map);
title('segmented Image:pp1')

figure;plot(tsmld,'*-b')
xlabel('Iteration');ylabel('Energy');
title('K-means energy minimization-pp1');

Answer 1

[Image Analyst]

I don't even see a call to kmeans() in your code. Adapt my attached example to however many spectral bands you want to use (it should be obvious how to do it).

[mohammed abdul wadood]

Thank you, actually I am just a beginner in matlab programming and with your help I am become better in matlab. thank you again.

[Image Analyst]

If you need help, attach all the spectral bands as a 3-D image, or as several separate grayscale images, in a .mat file. And tell me how many clusters you want to force there to be.

[mohammed abdul wadood]

Thank you sir, i have 7 sets of spectral image and i want k-means algorithm to segment any image in any set to 3 clusters just like cameramen example. I am sorry but i don't know how to attached all image in m.file, I can attached them from attached bottom above or send it on your email, and i well by grateful for helping my.

[Image Analyst]

You attached a demo gray scale image (covering all parts of the visible spectrum), not a spectral image. Please take your 7 spectral gray scale images (one monochrome image for each spectral wavelength band) and zip them together and attach.

Note, 7 spectral bands is NOT the same as one visible light image that you cluster into 7 gray level ranges (which is what you attached).

[mohammed abdul wadood]

yes sir what i attached before it's one of 7 gray scale image, i have image for two scenes, in the first scene i have 7 thermal band image and 14 blue band image and 7 gray scale image, in the second scene i have 7 visible image and 14 blue band image, that's all images i have it, i well zip all of them together and attached.

[Image Analyst]

Try this:

clc;    % Clear the command window.
close all;  % Close all figures (except those of imtool.)
clear;  % Erase all existing variables. Or clearvars if you want.
workspace;  % Make sure the workspace panel is showing.
format long g;
format compact;
fontSize = 25;

%===============================================================================
% Read in a demo image.
% Specify the folder where the files live.
myFolder = pwd;
% Check to make sure that folder actually exists.  Warn user if it doesn't.
if ~isdir(myFolder)
  errorMessage = sprintf('Error: The following folder does not exist:\n%s', myFolder);
  uiwait(warndlg(errorMessage));
  return;
end
% Get a list of all files in the folder with the desired file name pattern.
filePattern = fullfile(myFolder, 'pp*.bmp'); % Change to whatever pattern you need.
theFiles = dir(filePattern);
numberOfFiles = length(theFiles);
plotRows = ceil(sqrt(numberOfFiles));
for k = 1 : length(theFiles)
  baseFileName = theFiles(k).name;
  fullFileName = fullfile(myFolder, baseFileName);
  fprintf(1, 'Now reading %s\n', fullFileName);
  % Load image array with imread()
  spectralImage = imread(fullFileName);
  % It is supposed to be monochrome because it's a spectral image. Get the dimensions of the image.
  % numberOfColorChannels should be = 1 for a gray scale image, and 3 for an RGB color image.
  [rows, columns, numberOfColorChannels] = size(spectralImage);
  if numberOfColorChannels > 1
    % It's not really gray scale like we expected - it's color.
    % Use weighted sum of ALL channels to create a gray scale image.
    spectralImage = rgb2gray(spectralImage);
    % ALTERNATE METHOD: Convert it to gray scale by taking only the green channel,
    % which in a typical snapshot will be the least noisy channel.
    % grayImage = grayImage(:, :, 2); % Take green channel.
  end
  subplot(plotRows, plotRows, k);
  imshow(spectralImage);  % Display image.
  drawnow; % Force display to update immediately.
  title(baseFileName);

    % Load into columns of an array for kmeans()
    % Get the data for doing kmeans.  We will have one column for each color channel.
    if k == 1
      data = double(spectralImage(:));
    else
      data = [data, double(spectralImage(:))];
    end
  end
  % Enlarge figure to half screen.
  set(gcf, 'Units', 'Normalized', 'OuterPosition', [0, 0.04, 0.5, 0.96]);
  drawnow;
  hp = impixelinfo();

% %----------------------------------------------------------------------------------------
% % Ask user how many color classes they want.
% defaultValue = 6;
% titleBar = 'Enter an integer value';
% userPrompt = 'Enter the number of color classes to find ';
% caUserInput = inputdlg(userPrompt, titleBar, 1, {num2str(defaultValue)});
% if isempty(caUserInput),return,end % Bail out if they clicked Cancel.
% % Round to nearest integer in case they entered a floating point number.
% numberOfClasses = round(str2double(cell2mat(caUserInput)));
% % Check for a valid integer.
% if isnan(numberOfClasses) || numberOfClasses < 2 || numberOfClasses > 6
%   % They didn't enter a number.
%   % They clicked Cancel, or entered a character, symbols, or something else not allowed.
%   numberOfClasses = defaultValue;
%   message = sprintf('I said it had to be an integer.\nTry replacing the user.\nI will use %d and continue.', numberOfClasses);
%   uiwait(warndlg(message));
% end
numberOfClasses = 6; % Assume 6 for demo.

%----------------------------------------------------------------------------------------
%  KMEANS CLASSIFICATION RIGHT HERE!!!
% Each row of data represents one pixel.  Each column of data represents one image.
% Have kmeans decide which cluster each pixel belongs to.
indexes = kmeans(data, numberOfClasses);

%----------------------------------------------------------------------------------------
% Let's convert what class index the pixel is into images for each class index.
% Assume 6 clusters.
class1 = reshape(indexes == 1, rows, columns);
class2 = reshape(indexes == 2, rows, columns);
class3 = reshape(indexes == 3, rows, columns);
class4 = reshape(indexes == 4, rows, columns);
class5 = reshape(indexes == 5, rows, columns);
class6 = reshape(indexes == 6, rows, columns);
% Let's put these into a 3-D array for later to make it easy to display them all with a loop.
allClasses = cat(3, class1, class2, class3, class4, class5, class6);
allClasses = allClasses(:, :, 1:numberOfClasses); % Crop off just what we need.
% OK!  WE'RE ALL DONE!.  Nothing left now but to display our classification images.
figure;
plotRows = ceil(sqrt(size(allClasses, 3)));
% Display the classes, both binary and masking the original.
% Also make an indexes image so we can display each class in a unique color.
indexedImageK = zeros(rows, columns, 'uint8'); % Initialize another indexed image.
for c = 1 : numberOfClasses
  % Display binary image of what pixels have this class ID number.
  subplot(plotRows, plotRows, c);
  thisClass = allClasses(:, :, c);
  imshow(thisClass);
  caption = sprintf('Image of\nClass %d Indexes', c);
  title(caption, 'FontSize', fontSize);
end
% Enlarge figure to full screen.
set(gcf, 'Units', 'Normalized', 'OuterPosition', [0.5, 0.04, 0.5, 0.96]);

message = sprintf('Done with Classification');
helpdlg(message);

[uppu venkatesh]

can any one explain that the each image class indexes refers for , with name ? please can any one explain.

[ehsan goftari]

Hi Image analyst,

I want to cluster 3 images of band 2, NDVI and NDWI of landsat into four groups. can I use the code that you have commented here? will it work for me?

[Image Analyst]

Yes.

[Federico Alvarez]

Unrelated / related - thank you all (specially Image Analyst). I fell into a similar issue, and was suspecting that the way I was arranging my data was the problem. Looking at your code samples helped me confirm the problem.

Answer 2

[Jaswinder Singh]

can anybody provide me the code to find coefficient of a hysteresis loop whose coordinates are given to me. data is attached