function [Iseg,sep] = fun(I)
error(nargchk(1,2,nargin))
%% Checking n (number of classes)
if nargin==1
n = 2;
elseif n==1;
Iseg = NaN(size(I));
sep = 0;
return
elseif n~=abs(round(n)) || n==0
error('n must be a strictly positive integer!')
elseif n>255
n = 255;
warning(' n value has been changed to 255.')
end
I = single(I);
unI = sort(unique(I));
nbins = min(length(unI),256);
if nbins==n
Iseg = ones(size(I));
graycol = linspace(0,1,n);
for i = 1:n-1, Iseg(I==unI(i)) = graycol(i); end
sep = 1;
return
elseif nbins<n
Iseg = NaN(size(I));
sep = 0;
return
elseif nbins<256
[histo,pixval] = hist(I(:),unI);
else
[histo,pixval] = hist(I(:),256);
end
P = histo/sum(histo);
clear unI
w = cumsum(P);
mu = cumsum((1:nbins).*P);
if n==2
sigma2B =...
(mu(end)*w(2:end-1)-mu(2:end-1)).^2./w(2:end-1)./(1-w(2:end-1));
[maxsig,k] = max(sigma2B);
% segmented image
Iseg = ones(size(I),'uint8')*255;
Iseg(I<=pixval(k+1)) = 0;
% separability criterion
sep = maxsig/sum(((1:nbins)-mu(end)).^2.*P);
elseif n==3
w0 = w;
w2 = fliplr(cumsum(fliplr(P)));
[w0,w2] = ndgrid(w0,w2);
mu0 = mu./w;
mu2 = fliplr(cumsum(fliplr((1:nbins).*P))./cumsum(fliplr(P)));
[mu0,mu2] = ndgrid(mu0,mu2);
w1 = 1-w0-w2;
w1(w1<=0) = NaN;
sigma2B =...
w0.*(mu0-mu(end)).^2 + w2.*(mu2-mu(end)).^2 +...
(w0.*(mu0-mu(end)) + w2.*(mu2-mu(end))).^2./w1;
sigma2B(isnan(sigma2B)) = 0; % zeroing if k1 >= k2
[maxsig,k] = max(sigma2B(:));
[k1,k2] = ind2sub([nbins nbins],k);
% segmented image
Iseg = ones(size(I),'uint8')*255;
Iseg(I<=pixval(k1)) = 0;
Iseg(I>pixval(k1) & I<=pixval(k2)) = round(pixval(k2)/max(pixval)*255);
% separability criterion
sep = maxsig/sum(((1:nbins)-mu(end)).^2.*P);
else
k0 = linspace(0,1,n+1); k0 = k0(2:n);
[k,y] = fminsearch(@sig_func,k0);
k = round(k*(nbins-1)+1);
% segmented image
Iseg = ones(size(I),'uint8')*255;
Iseg(I<=pixval(k(1))) = 0;
maxpixval = max(pixval);
for i = 1:n-2
Iseg(I>pixval(k(i)) & I<=pixval(k(i+1))) =...
round(pixval(k(i+1))/maxpixval*255);
end
% separability criterion
sep = 1-y;
end
function y = sig_func(k)
muT = sum((1:nbins).*P);
sigma2T = sum(((1:nbins)-muT).^2.*P);
k = round(k*(nbins-1)+1);
k = sort(k);
if any(k<1 | k>nbins), y = 1; return, end
k = [0 k nbins];
sigma2B = 0;
for j = 1:n
wj = sum(P(k(j)+1:k(j+1)));
if wj==0, y = 1; return, end
muj = sum((k(j)+1:k(j+1)).*P(k(j)+1:k(j+1)))/wj;
sigma2B = sigma2B + wj*(muj-muT)^2;
end
y = 1-sigma2B/sigma2T; % within the range [0 1]
end
end
