clear, close all
B=imread('she.jpg');                                          %讀入原始jpg格式圖像
figure(1);
imshow(B),title('原始jpg格式圖像');
I1=rgb2gray(B);                                          %將原圖像轉化為灰度圖象
figure(2);
imshow(I1),title('灰度格式圖像');
BW1 = edge(I1,'sobel');
BW2 = edge(I1,'canny');
figure(6),imshow(BW1),title('邊緣檢測1');                                 %邊緣檢測
figure(5), imshow(BW2),title('邊緣檢測2');
[I1,map1]=gray2ind(I1,255);                              %將灰度圖像轉化為索引圖像
I1=double(I1);                                     %將unit8數組轉化為double型數組
 Z=I1                                             %將double型數組I1轉存到Z中
figure(3)                                              %畫出未進行分割的原始圖像
image(Z),title('未進行分割的原始圖像');colormap(map1);         
NIND=40;                                         %個體數目(Number of individuals)
MAXGEN=50;                        %最大遺傳代數(Maximum number of generations)
PRECI=8;                                  %變量的二進制位數(Precision of variables)
GGAP=0.9;                                                       %代溝(Generation gap)
FieldD=[8;1;256;1;0;1;1];                        %建立區域描述器(Build field descriptor)
Chrom=crtbp(NIND,PRECI);                                         %創建初始種群
gen=0;    
phen=bs2rv(Chrom,FieldD);                                    %初始種群十進制轉換
ObjV=target(Z,phen);                                           %計算種群適應度值
while gen<MAXGEN                                                %代溝(Generation gap)
    FitnV=ranking(-ObjV);                        %分配適應度值(Assign fitness values)
    SelCh=select('sus',Chrom,FitnV,GGAP);                                    %選擇
    SelCh=recombin('xovsp',SelCh,0.7);                                        %重組
    SelCh=mut(SelCh);                                                     %變異
    phenSel=bs2rv(SelCh,FieldD);                                  %子代十進制轉換
    ObjVSel=target(Z,phenSel);
    [Chrom ObjV]=reins(Chrom,SelCh,1,1,ObjV,ObjVSel);                      %重插入
    gen=gen+1;
end
[Y, I]=max(ObjV);
M=bs2rv(Chrom(I,:),FieldD);                                             %估計域值
[m, n]=size(Z);
for i=1:m
    for j=1:n
        if Z(i,j)>M                                      %灰度值大于域值時是白色
            Z(i,j)=256;
        end
    end
end
figure(4)                                                    %畫出分割后目標圖像
image(Z),title('分割后圖像');colormap(map1);

target求適應度函數代碼：
function f=target(T,M)                     %適應度函數,T為待處理圖像,M為域值序列
[U, V]=size(T);
W=length(M);
f=zeros(W,1);
for k=1:W
    I=0;s1=0;J=0;s2=0;                     %統計目標圖像和背景圖像的像素數及像素之和
    for i=1:U
        for j=1:V
            if T(i,j)<=M(k)
                s1=s1+T(i,j);I=I+1;
            end
            if T(i,j)>M(k)
                s2=s2+T(i,j);J=J+1;
            end
        end
    end
    if I==0,  p1=0;  else p1=s1/I; end
    if J==0,  p2=0;  else p2=s2/J; end
    f(k)=I*J*(p1-p2)*(p1-p2)/(256*256);
end;