%用otsu法對含噪圖象進行分割
C=imread('blwhite.jpg');
J=imnoise(C,'gaussian',0,400/(255*255));
subplot(1,3,1);imshow(C);title('原始灰度圖像');
subplot(1,3,2);imshow(J);title('加高斯噪聲的灰度圖像')
count=imhist(J);%直方圖統計
[r,t]=size(J); %圖象矩陣大小
N=r*t;%圖象象素個數
L=256;%指定圖像灰度級為256級
count=count/N;%各級灰度出現的概率
for i=1: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是灰度的起始值，Q是灰度的個數
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);%第tp個數最大，最大值為y
th=tp+p;%求出最大方差對應的灰度級
%if th<=110
 %   th=tp+p;
 
 %else th=110;%根據具體情況適當修正門限
 %end
Y1=zeros(r,t);
for i=1:r
    for j=1:t
        X1(i,j)=double(J(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
%以上代碼實現分割　
subplot(1,3,3);imshow(Y1);title('ostu灰度門限分割的圖像')