MATLAB实现图像分割otsuf 源程序代码_matlab

function main() %处理cameraman.tif I=imread('cameraman.tif'); %读入图像 Iseg=otsuf(I); %用otsuf程序分割 figure subplot(121) imshow(I) %绘制原图 subplot(122) imshow(Iseg) %绘制分割后的图 %处理circuit.tif I=imread('circuit.tif'); Iseg=otsuf(I); figure subplot(121) imshow(I) subplot(122) imshow(Iseg) %处理coins.png I=imread('coins.png'); Iseg=otsuf(I); figure subplot(121) imshow(I) subplot(122) imshow(Iseg) %处理lena.jpg I=imread('lena.jpg'); Iseg=otsuf(I); figure subplot(121) imshow(I) subplot(122) imshow(Iseg) %处理liftingbody.png I=imread('liftingbody.png'); Iseg=otsuf(I); figure subplot(121) imshow(I) subplot(122) imshow(Iseg) %处理rice.png I=imread('rice.png'); Iseg=otsuf(I); figure subplot(121) imshow(I) subplot(122) imshow(Iseg) %处理westconcordorthophoto.png I=imread('westconcordorthophoto.png'); Iseg=otsuf(I); figure subplot(121) imshow(I) subplot(122) imshow(Iseg) msgbox('MATLAB编程答疑，请加QQ: 1530497909','MATLAB答疑','help') web http://url.cn/TKcdXk -browser end function [Iseg,sep] = otsuf(I,n) %OTSU Gray-level image segmentation using Otsu's method. % Iseg = OTSU(I,n) computes a segmented image (Iseg) containing n classes % by means of Otsu's n-thresholding method (Otsu N, A Threshold Selection % Method from Gray-Level Histograms, IEEE Trans. Syst. Man Cybern. % 9:62-66;1979). Thresholds are computed to maximize a separability % criterion of the resultant classes in gray levels. % % Iseg = OTSU(I) uses n=2 (default value). The output Iseg is an unsigned % 8-bit integer image (i.e. of class uint8). % % [Iseg,sep] = OTSU(I,n) also returns the value (sep) of the separability % criterion within the range [0 1]. Zero is obtained only with images % having less than n gray level, whereas one (optimal value) is obtained % only with n-valued images. % % Notes: % ----- % It should be noticed that the thresholds generally become less credible % as the number of classes (n) to be separated increases (see Otsu's % paper for more details). % % Be careful! The OTSU function works with I of any size. An RGB image I % will be thus considered as a gray-level 3D-array! % % Example: % ------- % load clown % subplot(221) % X = ind2gray(X,map); % imshow(X) % title('Original','FontWeight','bold') % for n = 2:4 % Iseg = otsu(X,n); % subplot(2,2,n) % imshow(Iseg) % title(['n = ' int2str(n)],'FontWeight','bold') % end % % -- Damien Garcia -- 2007/08, revised 2009/03 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('MATLAB:otsu:WrongNValue',... 'n must be a strictly positive integer!') elseif n>255 n = 255; warning('MATLAB:otsu:TooHighN',... 'n is too high. n value has been changed to 255.') end %% Probability distribution 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 %% Zeroth- and first-order cumulative moments w = cumsum(P); mu = cumsum((1:nbins).*P); %% Maximal sigmaB^2 and Segmented image 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 to be minimized if n>=4 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