f.zip_zip

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