Why this program for automatical segmentation doesn´t work?

Dear all,

I have this code for automatical segmentation, but I don´t have output, MATLAB is only "busy".

function [mu,mask]=kmeans(ima,k)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%   kmeans image segmentation
%
%   Input:
%          ima: grey color image
%          k: Number of classes
%   Output:
%          mu: vector of class means 
%          mask: clasification image mask
%
%   Author: Jose Vicente Manjon Herrera
%    Email: [e-mail]
%     Date: 27-08-2005
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% check image
k=3;
ima = 'thorax-mdl.jpg';
ima=double(ima);
copy=ima;         % make a copy
ima=ima(:);       % vectorize ima
mi=min(ima);      % deal with negative 
ima=ima-mi+1;     % and zero values
s=length(ima);

% create image histogram
ima=round(ima); 
m=max(ima)+1;
h=zeros(1,m);
hc=zeros(1,m);
  

for i=1:s
  if(ima(i)>0) h(ima(i))=h(ima(i))+1;end;
end
ind=find(h);
hl=length(ind);

% initiate centroids
% k=str2double(k);
mu=(1:k)*m/(k+1);

% start process

while(true)

    oldmu=mu;
    % current classification  

    for i=1:hl
        c=abs(ind(i)-mu);
        cc=find(c==min(c));
        hc(ind(i))=cc(1);
    end

    %recalculation of means  

    for i=1:k, 
        a=find(hc==i);
        mu(i)=sum(a.*h(a))/sum(h(a));
    end

    if(mu==oldmu) break;end;

end

% calculate mask
s=size(copy);
mask=zeros(s);
for i=1:s(1),
for j=1:s(2),
  c=abs(copy(i,j)-mu);
  a=find(c==min(c));  
  mask(i,j)=a(1);
end
end

mu=mu+mi-1; 
imshow(mask,[])% recover real range

I don´t know, why I don´t have any output.I attach original image. Can you advise me?

Thank you for your answers.

Answer 1

[Jan]

Then obviously the condition if mu==oldmu is never true. This can happen even in the case of convergence if the result oscillate e.g. in the last sigfnificant digit.

You can monitor this easily:

count = 0;
while (true)
  ...
  count = count + 1;
  if count > 1000
    sprintf('%d iterations reached: mu = %.16g, oldmu = %.16g\n', count, mu, oldmu);
  end

    if (mu==oldmu) break; end
  end

Most likely a tolerance helps:

    if abs(mu - oldmu) < 1e-6
      break;
    end

[Veronika]

I transfered my code like this:

function [mu,mask]=kmeans(ima,k)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%   kmeans image segmentation
%
%   Input:
%          ima: grey color image
%          k: Number of classes
%   Output:
%          mu: vector of class means 
%          mask: clasification image mask
%
%   Author: Jose Vicente Manjon Herrera
%    Email: [e-mail]
%     Date: 27-08-2005
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% check image
k=3;
ima = 'thorax-mdl.jpg';
ima=double(ima);
copy=ima;         % make a copy
ima=ima(:);       % vectorize ima
mi=min(ima);      % deal with negative 
ima=ima-mi+1;     % and zero values
s=length(ima);

% create image histogram
ima=round(ima); 
m=max(ima)+1;
h=zeros(1,m);
hc=zeros(1,m);
  

for i=1:s
  if(ima(i)>0) h(ima(i))=h(ima(i))+1;end;
end
ind=find(h);
hl=length(ind);

% initiate centroids
% k=str2double(k);
mu=(1:k)*m/(k+1);

% start process
count = 0;%simon
while(true)

    oldmu=mu;
    % current classification  

    for i=1:hl
        c=abs(ind(i)-mu);
        cc=find(c==min(c));
        hc(ind(i))=cc(1);
    end

    %recalculation of means  

    for i=1:k, 
        a=find(hc==i);
        mu(i)=sum(a.*h(a))/sum(h(a));
    end

count = count + 1;
  if count > 1000
    sprintf('%d iterations reached: mu = %.16g, oldmu = %.16g\n', count, mu, oldmu);
  end
    if (mu==oldmu) break; end
  end
%   if(mu==oldmu) break;end;

end

% calculate mask
s=size(copy);
mask=zeros(s);
for i=1:s(1),
for j=1:s(2),
  c=abs(copy(i,j)-mu);
  a=find(c==min(c));  
  mask(i,j)=a(1);
end
end

mu=mu+mi-1; 
imshow(mask,[])% recover real range

But I have this error:

Error: File: kmeans.m Line: 79 Column: 1
This statement is not inside any function.
 (It follows the END that terminates the definition of the function "kmeans".)

I don´t know, where is mistake. Thank you for your answer.

[Jan]

Your code still contains "oldmu==mu". As explained already, tiny rounding errors will let this fail.

Use the auto-indentation to find the misplaced "end":

while(true)
  ...
  if (mu==oldmu) break; end
end  % ***
%   if(mu==oldmu) break;end;
end

Now the function kmeans() is finished. The next line

% calculate mask
s=size(copy);

Appears outside of any function, which is an error. I assume the end marked with * is wrong.

[Veronika]

Ok, so I removed this end with *, but still is Matlab only busy and no output...

function [mu,mask]=kmeans(ima,k)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%   kmeans image segmentation
%
%   Input:
%          ima: grey color image
%          k: Number of classes
%   Output:
%          mu: vector of class means 
%          mask: clasification image mask
%
%   Author: Jose Vicente Manjon Herrera
%    Email: [e-mail]
%     Date: 27-08-2005
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% check image
k=3;
ima = 'thorax-mdl.jpg';
ima=double(ima);
copy=ima;         % make a copy
ima=ima(:);       % vectorize ima
mi=min(ima);      % deal with negative 
ima=ima-mi+1;     % and zero values
s=length(ima);

% create image histogram
ima=round(ima); 
m=max(ima)+1;
h=zeros(1,m);
hc=zeros(1,m);
  

for i=1:s
  if(ima(i)>0) h(ima(i))=h(ima(i))+1;end;
end
ind=find(h);
hl=length(ind);

% initiate centroids
% k=str2double(k);
mu=(1:k)*m/(k+1);

% start process
count = 0;%simon
while(true)

    oldmu=mu;
    % current classification  

    for i=1:hl
        c=abs(ind(i)-mu);
        cc=find(c==min(c));
        hc(ind(i))=cc(1);
    end

    %recalculation of means  

    for i=1:k, 
        a=find(hc==i);
        mu(i)=sum(a.*h(a))/sum(h(a));
    end

count = count + 1;
  if count > 1000
    sprintf('%d iterations reached: mu = %.16g, oldmu = %.16g\n', count, mu, oldmu);
  end
    if (mu==oldmu) break; end

     if(mu==oldmu) break;end;

end

% calculate mask
s=size(copy);
mask=zeros(s);
for i=1:s(1),
for j=1:s(2),
  c=abs(copy(i,j)-mu);
  a=find(c==min(c));  
  mask(i,j)=a(1);
end
end

mu=mu+mi-1; 
imshow(mask,[])% recover real range

[Jan]

Sorry, Veronika, I'm telling you the third time now, that "mu==oldmu" cannot be expected to happen. The limited precision can lead to oscillating rounding effects even if te algorithm has converged. Inserting this command twice does not help also.

Please insert a tolerance:

if abs(mu - oldmu) < 1e-6  % Or whatever you accept as converged
  break;
end

[Veronika]

So I just removed rows, where mu==oldmu is and replaced this:

if abs(mu - oldmu) < 1e-6  % Or whatever you accept as converged
  break;
end

I´m so sorry, but I don´t understand, what you mean.