%利用最佳直方圖熵法（KSW熵法）實現灰度圖像閾值分割
%暫未涉及遺傳算法

filename=input('輸入預進行閾值分割的灰度圖像的文件名及路徑：','s');

disp('指定閾值搜索范圍（0～1）：');
lowvalue=input('下限值：');
highvalue=input('上限值：');

I=imread(filename);

info=imfinfo(filename);
maxgrayvalue=info.MaxSampleValue;
mingrayvalue=info.MinSampleValue;
graynum=maxgrayvalue-mingrayvalue+1;

hist=imhist(I);
total=0;
for i=mingrayvalue:maxgrayvalue
    total=total+hist(i+1);
end
hist1=hist/total;

HT=0;
for i=mingrayvalue:maxgrayvalue
    if hist1(i+1)==0
        temp=0;
    else
        temp=-hist1(i+1)*log(hist1(i+1));
    end
    HT=HT+temp;
end


%
lowvalue1=round(lowvalue*graynum);
highvalue1=round(highvalue*graynum);

for t=lowvalue1:highvalue1
    
    Pt=0;
    for i=1:t
        Pt=Pt+hist1(i);
    end
    
    Ht=0;
    for i=1:t
        if hist1(i)==0
           temp=0;
        else
           temp=-hist1(i)*log(hist1(i));
        end
        Ht=Ht+temp;
    end
    
    if Pt==0
        
    H(t)=log(Pt*(1-Pt))+Ht/Pt+(HT-Ht)/(1-Pt);
    
end

H(mingrayvalue+1:lowvalue1-1)=0;
H(highvalue1+1:maxgrayvalue+1)=0;

max_shannon=max(H);
t_opt=find(H==max_shannon);

threshold_opt=t_opt/graynum;

I1=im2bw(I,threshold_opt);

disp('灰度圖像閾值分割的效果如圖所示：');
disp('源圖為：Fifure No.1');
disp('最佳直方圖熵法閾值分割后的圖像為：Fifure No.2');
disp('Otsu法閾值分割后的圖像為：Fifure No.3');

figure(1);
imshow(I);
title('源圖');

figure(2);
imshow(I1);
title('最佳直方圖熵法閾值分割后的圖像');

level=graythresh(I);
I2=im2bw(I,level);
figure(3);
imshow(I2);
title('Otsu法閾值分割后的圖像');