利用最佳直方圖熵法（KSW熵法）及傳統遺傳演算法實現灰度影象二閾值分割

matlab程式碼如下：
1、main.m(主函式)：

%%%利用最佳直方圖熵法（KSW熵法）及傳統遺傳演算法實現灰度影象二閾值分割
%%%主程式
%%  初始部分，讀取影象及計算相關資訊
 clear;
 close all;
 clc;
I=imread('D:\MATLAB\work\2.21.jpg');
figure
figure(1),imshow(I);
 I=rgb2gray(I);


% I=imread('Lenna.bmp');

hist=imhist(I);     %顯示影象資料柱狀圖
 

total=0;
for i=0:255
    total=total+hist(i+1);
end
hist1=hist/total;   %求每點的歸一化//求畫素為i的概率Pi


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% 程式主幹部分

    %種群隨機初始化，種群數取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
        if X0(i,1)>X0(i,2)
            temp=X0(i,1);
            X0(i,2)=temp;
            X0(i,1)=temp;
        else
        end
        adapt_value0(i)=ksw_2(X0(i,1),X0(i,2),0,255,hist1);
    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_2d(X1,adapt_value1);

        s_code10=dec2bin(s1(:,1),8);
        s_code11=dec2bin(s1(:,2),8);

        [c10,c11]=cross_2d(s_code10,s_code11);

        [v10,v11]=mutation_2d(c10,c11);

        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_2(X2(i,1),X2(i,2),0,255,hist1);
        end

        adapt_average2=mean(adapt_value2);

        if abs(adapt_average2-adapt_average1)<=0.03
            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),:);

    t1=clock;
    search_time=etime(t1,t0);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%  閾值分割及顯示部分

I_temp1=I;
[height,width]=size(I_temp1)
for i=1:height
    for j=1:width
         if I_temp1(i,j)<opt(1);
                 I_temp1(i,j)=0;
        else if I_temp1(i,j)>opt(2);
                 I_temp1(i,j)=255;
            else I_temp1(i,j)=180;
            end
        end
    end
end
I1= I_temp1;
disp('灰度影象閾值分割的效果如圖所示：');
disp('源圖為：Fifure No.1');
disp('最佳直方圖熵法及傳統遺傳演算法閾二值分割後的影象為：Fifure No.2');

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

figure(3);
imshow(I1);
title('最佳直方圖熵法及傳統遺傳演算法閾二值分割後的影象');


disp('最佳直方圖熵法及傳統遺傳演算法二閾值為(s,t)：');
disp(opt(1));
disp(opt(2));

disp('最佳直方圖熵法及傳統遺傳演算法二閾值搜尋所用時間（s）：');
disp(search_time);

%%  程式結束

2、子函式

function s1=select_2d(X1,adapt_value1)

    %選擇運算元

    population=20;

    total_adapt_value1=0;
    for i=1:population
        total_adapt_value1=total_adapt_value1+adapt_value1(i);
    end
    adapt_value1_new=adapt_value1/total_adapt_value1;

    r=rand(1,population);

    for i=1:population
        temp=0;
        for j=1:population
            temp=temp+adapt_value1_new(j);
            if temp>=r(i)
                s1(i,:)=X1(j,:);
                break;
            end
        end
    end

function [c10,c11]=cross_2d(s_code10,s_code11)

   %交叉運算元

   pc=0.8;       %交叉概率取0.6
   population=20;

   %(1,2)/(3,4)/(5,6)進行交叉運算，(7,8)/(9,10)複製

   ww0=s_code10;
   ww1=s_code11;

   for i=1:(pc*population/2)
       r0=abs(round(rand(1)*10)-3);
       r1=abs(round(rand(1)*10)-3);
       for j=(r0+1):8
           temp0=ww0(2*i-1,j);
           ww0(2*i-1,j)=ww0(2*i,j);
           ww0(2*i,j)=temp0; 
       end
       for j=(r1+1):8
           temp1=ww1(2*i-1,j);
           ww1(2*i-1,j)=ww1(2*i,j);
           ww1(2*i,j)=temp1; 
       end
   end

   c10=ww0;
   c11=ww1;

function [v10,v11]=mutation_2d(c10,c11)

    %變異運算元

    format long;

    population=20;

    pm=0.03;

    for i=1:population
        for j=1:8
            r0=rand(1);
            r1=rand(1);
            if r0>pm
                temp0(i,j)=c10(i,j);
            else
                tt=not(str2num(c10(i,j)));
                temp0(i,j)=num2str(tt);
            end
            if r1>pm
                temp1(i,j)=c11(i,j);
            else
                tt=not(str2num(c11(i,j)));
                temp1(i,j)=num2str(tt);
            end
        end
    end

    v10=temp0;
    v11=temp1;

function y=ksw_2(s,t,mingrayvalue,maxgrayvalue,hist1)


   %計算最佳直方圖熵（KSW熵）

    Ps=0;%初始化
    for i=mingrayvalue:s; %從0到s
        Ps=Ps+hist1(i+1);%求和
    end

    Pt=0;%初始化
    for i=s:t; %從s+1到t
        Pt=Pt+hist1(i+1);%求和
    end
     Pn=0;%初始化
    for i=t:maxgrayvalue; %從t+1到n
        Pn=Pn+hist1(i+1);%求和
    end
    Hs=0;
    for i=mingrayvalue:s
        if hist1(i+1)==0%直方圖值為零者賦零
           temp=0;
        else
           temp=hist1(i+1)*log(1/hist1(i+1));%
        end
        Hs=Hs+temp;%
    end    
     Ht=0;
    for i=s:t
        if hist1(i+1)==0%直方圖值為零者賦零
           temp=0;
        else
           temp=hist1(i+1)*log(1/hist1(i+1));%
        end
        Ht=Ht+temp;%
    end  
     Hn=0;
    for i=t:maxgrayvalue
        if hist1(i+1)==0%直方圖值為零者賦零
           temp=0;
        else
           temp=hist1(i+1)*log(1/hist1(i+1));%
        end
        Hn=Hn+temp;%
    end

    if Ps==0 || Ps==1||Pt==0 || Pt==1||Pn==0 || Pn==1         %   or(Pt==0,Pt==1)
        temp1=0;
    else 
        temp1=log(Ps)+log(Pt)+log(Pn)+Hs/Ps+Ht/Pt+Hn/Pn;%影象總熵
    end

    if temp1 < 0
        H=0;
    else
        H=temp1;
    end


    y=H;

執行結果如下：

原圖：

灰度影象：

雙閾值分割影象：