%***************************************************************************************************************************************************************
%OTSU法对于具有双峰性质的灰度图像或是彩色图像的某一通道的分割效果很好,程序为了增加健壮性加了个可以根据实际情况确定的修正值th_set.
%**************************************************************************************************************************************************************
function Threshold=OTSU_Origin(image)
%a2=imread('color1.bmp');
%gray=rgb2gray(image);%原图像的灰度图
gray=uint8(image); %转换成图像的矩阵数值。
%figure;imshow(gray);
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;
st=0;
nd=254;
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和q分别是灰度的起始和结束值
u=0;
ua=0;
for i=1:q;
u=u+f(i)*(p+i-1);
ua(i)=u;
end
%计算图像的平均灰度值 %求均值向量sigma(i*p)
w=0.5;
for i=1:q;
w(i)=sum(f(1:i));
end
%计算出选择不同k的时候,A区域的概率
d=(u*w-ua).^2./(w.*(1-w));%求出不同k值时类间方差
%tp=max(d)
%[y,tp]=max(d)%求出最大方差对应的灰度级
Threshold=max(d);%求出最大方差
%Threshold=tp+p;
%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('灰度门限分割的图像');