基于人工蜂群的二维OTSU算法

已经过实验，可运行

人工蜂群算法部分

%人工蜂群算法部分
%experiment.m
%人工蜂群算法第一次应用
clear all                                     %preprocessing预处理
close all
clc
%图像预处理
%I=imread('E:\lena512.bmp');
%R=imread('E:\car.bmp');
%R=imread('E:\bird.jpg');
%I=imread('AT3_1m4_01.tif');
%I=imread('rice.png');
R=imread('peppers.png');
I=rgb2gray(R);
%二维直方图平面分布图
row=size(I,1);
column=size(I,2);
N0=row*column;
M=zeros(256,256);%区域平均灰度值图像
for i=1:row
    for j=1:column
        k=0;
        for i1=-1:1
            for i2=-1:1
                try %边缘区域处理
                    k=k+I(i+i1,j+i2)/9;
                catch
                    %k=k;
                end
            end
        end%邻域平均值计算
        if (I(i,j)==0||k==0)
        else
        M(I(i,j),k)=M(I(i,j),k)+1;%记录用以生成像素周边区域平均像素
        M1(i,j)=k;
        end
    end  %图像循环
end
 %ui的计算
    sum0=0;
    for s1=1:256
        for t1=1:256
            sum0=sum0+M(s1,t1)*s1;
        end
    end
    ui=sum0/N0;
    %uj的计算
    sum0=0;
    for s1=1:256
        for t1=1:256
            sum0=sum0+M(s1,t1)*t1;
        end
    end
    uj=sum0/N0;
    u=[ui,uj];
%人工蜂群算法部分
NP=50; %蜂群规模
FoodNumber=NP/2; %食物源数量
limit=100; %/*A food source which could not be improved through "limit" trials is abandoned by its employed bee*/蜂群规模*维数比较合适
maxCycle=600; %觅食周期数--循环次数
%参数定义
objfun='DDOTSU'; %cost function to be optimized    成本函数
D=4; %/*The number of parameters of the problem to be optimized   要优化参数数量*/
ub=ones(1,D)*255; %/*lower bounds of the parameters.   参数下界*/
lb=ones(1,D)*(1);%/*upper bound of the parameters.*/

runtime=1;%/*Algorithm can be run many times in order to see its robustness*/鲁棒性检测

GlobalMins=zeros(1,runtime);

for r=1:runtime
 
% /*All food sources are initialized */解空间生成
%/*Variables are initialized in the range [lb,ub]. If each parameter has different range, use arrays lb[j], ub[j] instead of lb and ub */

Range = repmat((ub-lb),[FoodNumber 1]);
Lower = repmat(lb, [FoodNumber 1]);
Foods = rand(FoodNumber,D) .* Range + Lower;
Foods=round(Foods);

ObjVal=feval(objfun,Foods,M,N0,u);
Fitness=calculateFitness(ObjVal);

%reset trial counters
trial=zeros(1,FoodNumber);

%/*The best food source is memorized*/
BestInd=find(ObjVal==min(ObjVal));
BestInd=BestInd(end);
GlobalMin=ObjVal(BestInd);
GlobalParams=Foods(BestInd,:);

iter=1;
while ((iter <= maxCycle)),

%%%%%%%%% EMPLOYED BEE PHASE 雇佣蜂%%%%%%%%%%%%%%%%%%%%%%%%
    for i=1:(FoodNumber)
        
        %/*The parameter to be changed is determined randomly*/
        Param2Change=fix(rand*D)+1;
        
        %/*A randomly chosen solution is used in producing a mutant solution of the solution i*/
        neighbour=fix(rand*(FoodNumber))+1;
       
        %/*Randomly selected solution must be different from the solution i*/        
            while(neighbour==i)
                neighbour=fix(rand*(FoodNumber))+1;
            end;
        
       sol=Foods(i,:);
       %  /*v_{ij}=x_{ij}+\phi_{ij}*(x_{kj}-x_{ij}) */
       sol(Param2Change)=Foods(i,Param2Change)+(Foods(i,Param2Change)-Foods(neighbour,Param2Change))*(rand-0.5)*2;
       sol=round(sol);
       %  /*if generated parameter value is out of boundaries, it is shifted onto the boundaries*/
        ind=find(sol<lb);
        sol(ind)=lb(ind);
        ind=find(sol>ub);
        sol(ind)=ub(ind);
        
        %evaluate new solution
        ObjValSol=feval(objfun,sol,M,N0,u);
        FitnessSol=calculateFitness(ObjValSol);
        
       % /*a greedy selection is applied between the current solution i and its mutant*/
       if (FitnessSol>Fitness(i)) %/*If the mutant solution is better than the current solution i, replace the solution with the mutant and reset the trial counter of solution i*/
            Foods(i,:)=sol;
            Fitness(i)=FitnessSol;
            ObjVal(i)=ObjValSol;
            trial(i)=0;
        else
            trial(i)=trial(i)+1; %/*if the solution i can not be improved, increase its trial counter*/
       end;
         
    end;

%%%%%%%%%%%%%%%%%%%%%%%% CalculateProbabilities %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%/* A food source is chosen with the probability which is proportioal to its quality*/
%/*Different schemes can be used to calculate the probability values*/
%/*For example prob(i)=fitness(i)/sum(fitness)*/
%/*or in a way used in the metot below prob(i)=a*fitness(i)/max(fitness)+b*/
%/*probability values are calculated by using fitness values and normalized by dividing maximum fitness value*/

prob=(0.9.*Fitness./max(Fitness))+0.1;
  
%%%%%%%%%%%%%%%%%%%%%%%% ONLOOKER BEE PHASE %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

i=1;
t=0;
while(t<FoodNumber)
    if(rand<prob(i))
        t=t+1;
        %/*The parameter to be changed is determined randomly*/
        Param2Change=fix(rand*D)+1;
        
        %/*A randomly chosen solution is used in producing a mutant solution of the solution i*/
        neighbour=fix(rand*(FoodNumber))+1;
       
        %/*Randomly selected solution must be different from the solution i*/        
            while(neighbour==i)
                neighbour=fix(rand*(FoodNumber))+1;
            end;
        
       sol=Foods(i,:);
       %  /*v_{ij}=x_{ij}+\phi_{ij}*(x_{kj}-x_{ij}) */
       sol(Param2Change)=Foods(i,Param2Change)+(Foods(i,Param2Change)-Foods(neighbour,Param2Change))*(rand-0.5)*2;
       sol=round(sol);        
       %  /*if generated parameter value is out of boundaries, it is shifted onto the boundaries*/
        ind=find(sol<lb);
        sol(ind)=lb(ind);
        ind=find(sol>ub);
        sol(ind)=ub(ind);
        
        %evaluate new solution

        ObjValSol=feval(objfun,sol,M,N0,u);
        FitnessSol=calculateFitness(ObjValSol);
       
       % /*a greedy selection is applied between the current solution i and its mutant*/
       if (FitnessSol>Fitness(i)) %/*If the mutant solution is better than the current solution i, replace the solution with the mutant and reset the trial counter of solution i*/
            Foods(i,:)=sol;
            Fitness(i)=FitnessSol;
            ObjVal(i)=ObjValSol;
            trial(i)=0;
        else
            trial(i)=trial(i)+1; %/*if the solution i can not be improved, increase its trial counter*/
       end;
    end;
    
    i=i+1;
    if (i==(FoodNumber)+1) 
        i=1;
    end;   
end; 
%/*The best food source is memorized*/
         ind=find(ObjVal==min(ObjVal));
         ind=ind(end);
         if (ObjVal(ind)<GlobalMin)
         GlobalMin=ObjVal(ind);
         GlobalParams=Foods(ind,:);
         end;
%%%%%%%%%%%% SCOUT BEE PHASE %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%/*determine the food sources whose trial counter exceeds the "limit" value. 
%In Basic ABC, only one scout is allowed to occur in each cycle*/

ind=find(trial==max(trial));
ind=ind(end);
if (trial(ind)>limit)
    Bas(ind)=0;
    sol=(ub-lb).*rand(1,D)+lb;
    sol=round(sol);
    ObjValSol=feval(objfun,sol,M,N0,u);
    FitnessSol=calculateFitness(ObjValSol);
    Foods(ind,:)=sol;
    Fitness(ind)=FitnessSol;
    ObjVal(ind)=ObjValSol;
end;

%fprintf('iter=%d ObjVal=%g\n',iter,GlobalMin);
iter=iter+1;

end % End of ABC

GlobalMins(r)=GlobalMin;
end; %end of runs
%fprintf('pos=%g\n',mean(pos,2));
%fprintf('mean=%g\n',mean(GlobalMins));
save all

%人工蜂群算法部分结束
BestValue=[(GlobalParams(1)+GlobalParams(3))/2 (GlobalParams(2)+GlobalParams(4))/2];
if BestValue(1)>BestValue(2)
    cach=BestValue(1);
    BestValue(1)=BestValue(2);
    BestValue(2)=cach;
end
%根据阈值进行图像分割
for i=1:row
    for j=1:column
        if(I(i,j)>=BestValue(1)&&I(i,j)<=BestValue(2))
            K(i,j)=1;
        else
            K(i,j)=0;
        end
    end
end
figure;
subplot(121);imshow(I);
subplot(122);imshow(K);

二维双阈值算法部分

%DDOTSU.m
function ObjVal = DDOTSU(Chrom,M,N0,u)
FoodNumber=size(Chrom,1);
ObjVal=zeros(1,FoodNumber);
for fn=1:FoodNumber
    s1=Chrom(fn,1);
    s2=Chrom(fn,2);
    t1=Chrom(fn,3);
    t2=Chrom(fn,4);
%%
    %w0的计算
    sum0=0;
    for s=s1:s2
        for t=t1:t2
            sum0=sum0+M(s,t);
        end
    end
    w0=sum0/N0;
    %w1的计算
    sum0=0;
    for s=1:s1
        for t=1:t2
            sum0=sum0+M(s,t);
        end
    end
    w1=sum0/N0;
    sum0=0;
    for s=s2:256
        for t=t2:256
            sum0=sum0+M(s,t);
        end
    end
    w2=sum0/N0;
%%
    %ui0的计算
    sum0=0;
    for s=s1:s2
        for t=t1:t2
            sum0=sum0+M(s,t)*s;
        end
    end
    ui0=sum0/N0;
    %uj0的计算
    sum0=0;
    for s=s1:s2
        for t=t1:t2
            sum0=sum0+M(s,t)*t;
        end
    end
    uj0=sum0/N0;
    u0=[ui0,uj0];
%%
    %ui1的计算
    sum0=0;
    for s=1:s1
        for t=1:t1
            sum0=sum0+M(s,t)*s;
        end
    end
    ui1=sum0/N0;
    sum0=0;
    for s=s2:256
        for t=t2:256
            sum0=sum0+M(s,t)*s;
        end
    end
    ui2=sum0/N0;
    %uj1的计算
    sum0=0;
    for s=1:s1
        for t=1:t1
            sum0=sum0+M(s,t)*t;
        end
    end
    uj1=sum0/N0;
    sum0=0;
    for s=s2:256
        for t=t2:256
            sum0=sum0+M(s,t)*t;
        end
    end
    uj2=sum0/N0;
    u1=[ui1,uj1];
    u2=[ui2,uj2];
%%
    trsb=w0*(u0-u)*(u0-u)'+w1*(u1-u)*(u1-u)'+w2*(u2-u)*(u2-u)';
    ObjVal(fn)=1/trsb;
end
end

评估函数

%calculateFitness.m
function fFitness=calculateFitness(fObjV)
fFitness=zeros(size(fObjV));
ind=find(fObjV>=0);
fFitness(ind)=1./(fObjV(ind)+1);
ind=find(fObjV<0);
fFitness(ind)=1+abs(fObjV(ind));

原文地址：https://www.cnblogs.com/zhoushuaiyi/p/14670703.html