%%%利用二維最佳直方圖熵法（KSW熵法）及窮舉法實現灰度圖像閾值分割
%%%主程序


%%初始部分，讀取圖像及計算相關信息

clear;
close all;
clc;

%format long;

I=imread('rice.tif');


windowsize=3;
I_temp=I;
for i=2:255
    for j=2:255
        I_temp(i,j)=round(mean2(I(i-1:i+1,j-1:j+1)));
    end
end
I_average=I_temp;


I_p=I;
I_average_p=I_average;
hist_2d(1:256,1:256)=zeros(256,256);
for i=1:256
    for j=1:256
        hist_2d(I_p(i,j),I_average_p(i,j))=hist_2d(I_p(i,j),I_average_p(i,j))+1;
    end
end

total=256*256;

hist_2d_1=hist_2d/total;


%%%%%%

Hst=0;
for i=0:255
    for j=0:255
        if hist_2d_1(i+1,j+1)==0
            temp=0;
        else
            temp=hist_2d_1(i+1,j+1)*log(1/hist_2d_1(i+1,j+1));
        end
        Hst=Hst+temp;
    end
end



%%程序主干部分
t0=clock;

    for s=0:255
        for t=0:255
            adapt_value(s+1,t+1)=ksw_2d(s,t,0,255,hist_2d_1,Hst);
        end
    end
        
    
    [max_value1,index1]=max(adapt_value);
    [max_value2,index2]=max(max_value1);
    t_opt=index2-1;
    s_opt=index1(index2)-1;
    
t1=clock;
search_time=etime(t1,t0);
    
%%閾值分割及顯示部分

threshold_opt=s_opt/255;

I1=im2bw(I,threshold_opt);

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

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

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


disp('二維最佳直方圖熵法及窮舉法閾值為(s,t)：');
disp(s_opt);
disp(t_opt);

disp('二維最佳直方圖熵法及窮舉法閾值搜索所用時間（s）：');
disp(search_time);

%%程序結束