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


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


clear;
close all;
clc;

%format long;

I=imread('Lenna_noise.bmp');


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=imadd(I,1);
I_average_p=imadd(I_average,1);
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



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%程序主干部分

    %種群隨機初始化，種群數取20，染色體二進制編碼取16位
    t0=clock;
    
    population=20;
  
    X00=round(rand(1,population)*255);
    X01=round(rand(1,population)*255);
    
    for i=1:population
        X0(i,:)=[X00(i) X01(i)];
    end
    
    
    for i=1:population
        adapt_value0(i)=ksw_2d(X0(i,1),X0(i,2),0,255,hist_2d_1,Hst);
    end
    
    adapt_average0=mean(adapt_value0);
    
    X1=X0;
    adapt_value1=adapt_value0;
    adapt_average1=adapt_average0;
    
    %循環搜索，搜索代數取100
    
    generation=100;
    
    for k=1:generation
        
        s1=select_2d1(X1,adapt_value1);
    
        s_code10=dec2bin(s1(:,1),8);
        s_code11=dec2bin(s1(:,2),8);
    
        [c10,c11]=cross_2d1(s_code10,s_code11,k);
    
        [v10,v11]=mutation_2d1(c10,c11,k);
    
        X20=(bin2dec(v10))';
        X21=(bin2dec(v11))';
        
        for i=1:population
            X2(i,:)=[X20(i) X21(i)];
        end
    
        for i=1:population
            adapt_value2(i)=ksw_2d(X2(i,1),X2(i,2),0,255,hist_2d_1,Hst);
        end
    
        adapt_average2=mean(adapt_value2);
    
        if abs(adapt_average2-adapt_average1)<=0.072
            break;
        else
            X1=X2;
            adapt_value1=adapt_value2;
            adapt_average1=adapt_average2; 
        end
    end
    
    max_value=max(adapt_value2);
    number=find(adapt_value2==max_value);
    opt=X2(number(1),:);
    s_opt=opt(1);
    t_opt=opt(2);

    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);

%%程序結束