***************************************************************************************************************************************************************

OTSU法對于具有雙峰性質的灰度圖像或是彩色圖像的某一通道的分割效果很好,程序為了增加健壯性加了個可以根據實際情況確定的修正值th_set.

**************************************************************************************************************************************************************


function y1=OTSU(image,th_set)

%a2=imread('color1.bmp');
gray=rgb2gray(image);%原圖像的灰度圖
low_high=stretchlim(gray);%增強圖像，似乎也不是一定需要
gray=imadjust(gray,low_high,[]);
% subplot(224);imshow(gray);title('after adjust');

count=imhist(gray);
[r,t]=size(gray);
n=r*t;
l=256;
count=count/n;%各級灰度出現的概率
for i=2:l
        if count(i)~=0
            st=i-1;
            break
        end
end
%以上循環語句實現尋找出現概率不為0的最小灰度值
for i=l:-1:1
        if count(i)~=0;
            nd=i-1;
            break
        end
end
%實現找出出現概率不為0的最大灰度值
f=count(st+1:nd+1);
p=st;q=nd-st;%p和分別是灰度的起始和結束值
u=0;
for i=1:q;
        u=u+f(i)*(p+i-1);
        ua(i)=u;
end
%計算圖像的平均灰度值
for i=1:q;
        w(i)=sum(f(1:i));
end
%計算出選擇不同k的時候，A區域的概率
d=(u*w-ua).^2./(w.*(1-w));%求出不同k值時類間方差
[y,tp]=max(d);%求出最大方差對應的灰度級
th=tp+p;
if th<th_set
        th=tp+p;
else
        th=th_set;           %根據具體情況適當修正門限
end
y1=zeros(r,t);
for i=1:r
        for j=1:t
            x1(i,j)=double(gray(i,j));
        end
end
for i=1:r
        for j=1:t
            if (x1(i,j)>th)
                y1(i,j)=x1(i,j);
            else
                y1(i,j)=0;
            end
        end
end
%上面一段代碼實現分割
% figure,imshow(y1);
% title('灰度門限分割的圖像');