frew.m

ostu图像分割阈值算法,对Ostu图像分割最优阈值进行优化处理

clear;

%灰度图象的直方图数据

fn='lena.bmp';
[I,map]=imread (fn);
figure;imshow(I);
J=I;
global G_min G_max
G_min =double(min(min(I)));    % the min gray value of the image
G_max =double(max(max(I)));    % the max gray value of the iamge

%计算灰度图象的直方图数据
L=256; %灰度级 
global p
p=zeros(L,1);%统计直方图结果数据
nk=zeros(L,1);
[row,col]=size(I);

s=row*col; %总像素个数
for i = 1:row
for j = 1:col
num = double(I(i,j))+1; %获取像素点灰度级
nk(num,1) = nk(num,1)+1; %统计nk 
end
end
%计算直方图概率估计
for i=1:L
p(i,1)=nk(i,1)/s; 
end

MaxCycle = 100;%循环次数
SltnDmsnNmbr = 1;%维数

SltnNmbr = 4;%竞选人个数

SltnMaxLmt = G_max;%最大值，最小值范围
SltnMinLmt = G_min;

Invst2SltnNmbr = 1;%？？

SltnLclInvstNmbr = 3;%正态分布个数
SltnLclInvstVar = 0.02;

SltnSpcInvstNmbr = 12;%均匀分布个数

SltnInvstNmbr = SltnNmbr * SltnLclInvstNmbr + SltnSpcInvstNmbr;%总的个数

for i = 1 : SltnNmbr
    Sltn(1, i, :) = GnrtSpcSltn(SltnDmsnNmbr, SltnMaxLmt, SltnMinLmt);%产生竞选人的位置
end
%CmptCycl=1;
for CmptCycle = 1 : MaxCycle%循环次数
   
    for i = 1 : SltnNmbr
        for j = 1 : SltnLclInvstNmbr
            Invst((i - 1) * SltnLclInvstNmbr + j, :) = GnrtLclInvst(Sltn(CmptCycle, i, :), SltnLclInvstVar, SltnMaxLmt, SltnMinLmt);%竞选人周围正态分布
        end
    end

    for i = 1 : SltnSpcInvstNmbr
        Invst(SltnNmbr * SltnLclInvstNmbr + i, :) = GnrtSpcSltn(SltnDmsnNmbr, SltnMaxLmt, SltnMinLmt);%全局均匀分布
    end

    for i = 1 : SltnInvstNmbr
        InvstValue(i) = ostu(Invst(i, :));%计算所有人的适应值
    end

    for i = 1 : SltnNmbr
        for j = 1 : SltnInvstNmbr
            TempSum = 0;
            for k = 1 : SltnDmsnNmbr
                TempSum = TempSum + abs(Invst(j, k) - Sltn(CmptCycle, i, k));
            end
            Dstnc(i, j) = TempSum;
        end         
    end

    for i = 1 : SltnNmbr
        sum(i) = 0;
        for j = 1 : SltnInvstNmbr
            sum(i) = sum(i) + InvstValue(i) / Dstnc(i,j);    
        end
    end

    for i = 1 : SltnNmbr
        for j = 1 : SltnInvstNmbr
            Power(i, j) = InvstValue(i) / Dstnc(i,j) / sum(i);    
        end
    end

    for i = 1 : SltnNmbr   
        for k = 1 : SltnDmsnNmbr
            ThisSum = 0;
            for j = 1 : SltnInvstNmbr
                ThisSum = ThisSum + Power(i,j) * Invst(j,k);   
            end
            Sltn(CmptCycle + 1, i, k)=ThisSum;
        end    
    end    
    
    for i = 1 : SltnNmbr
        SltnValue(i) = ostu(Sltn(CmptCycle + 1, i, :)); 
    end

    for k = 1 : Invst2SltnNmbr
        
        MinSltnValue = min(SltnValue);
        
        for i1 = 1 : SltnNmbr
            if SltnValue(i1) == MinSltnValue
                MinSltn = Sltn(CmptCycle+1,i1,:); 
                break;
            end    
        end

        MaxInvstValue = max(InvstValue);

        for i2 = 1 : SltnInvstNmbr
            if InvstValue(i2) == MaxInvstValue
                MaxInvst = Invst(i2,:);
                break;
            end
        end

        if MaxInvstValue > MinSltnValue
            Sltn(CmptCycle+1, i1, :) = MaxInvst;
            SltnValue(i1) = MaxInvstValue;
            Invst(i2, :) = MinSltn;
            InvstValue(i2) = MinSltnValue;
        end
    end

end

SltnWndwHigh = 1.2 * max(max(max(Sltn)));

for i = 1 : SltnDmsnNmbr
    figure; 
    plot(Sltn(:, :, i));
    axis([0, MaxCycle, SltnMinLmt, SltnMaxLmt]);
end

for i = 1 : MaxCycle
    for j = 1 : SltnNmbr
        x(i) = i;
        y(j, i) = ostu(Sltn(i,j,:));
    end
end

figure;
plot(x,y);
axis([0, MaxCycle, 0,10000]);
figure;
BW=im2BW(I,map,round(MaxInvst)/255);
imshow(BW);