基于万有引力搜索算法图像分割的MATLAB实现

处理效果：

1. 原图

2. 处理结果

3. 相关参数

种群规模：5
种群最大迭代次数：20
万有引力算法计算出的阈值：156.2703
关于万有引力算法的程序代码都来自http://blog.csdn.net/u013337691/article/details/52732631
以下为具体程序代码：

1. 图像处理相关程序

%% 清空环境变量
close all
clear
clc
format compact
%% 选择图片，并二值化
[fn,pn,fi]=uigetfile('*.jpg','选择图片');
I=imread([pn fn]); 
if ndims(I) == 3
I = rgb2gray(I);
end
% fxy = imhist(I, 256); %统计每个灰度值的个数
[counts,x] = imhist(I, 256) ;
figure;
subplot(2, 2, 1); 
imshow(I, []); title('原图')

%% GSA优化参数 
N=5; % 群体规模 Number of agents.
max_it=10; % 最大迭代次数 Maximum number of iterations (T).
ElitistCheck=0; % 如果ElitistCheck=1,则使用文献中的公式21；如果ElitistCheck=0，则用文献中的公式9.
Rpower=1;% 文献中公式7中的R的幂次数 power of 'R' in eq.7.
min_flag=0; % 取1求解极小值问题，取0求解极大值问题 1: minimization, 0: maximization.
objfun=@objfun_image; % 目标函数
[Fbest,Lbest,BestChart,MeanChart]=GSA_image(objfun,N,max_it,ElitistCheck,min_flag,Rpower,...
    counts,x);
Fbest;
Lbest

p=Lbest(1)/255;


% Fbest: 最优目标值 Best result. 
% Lbest: 最优解 Best solution. The location of Fbest in search space.
% BestChart: 最优解变化趋势 The best so far Chart over iterations. 
% MeanChart: 平均适应度函数值变化趋势 The average fitnesses Chart over iterations.

%subplot(2, 2, 2);
%plot(fxy); %画出灰度直方图
%title('直方图')
% p1 = {'Input Num:'}; 
% p2 = {'180'};   %手动输入阈值
% p3 = inputdlg(p1,'Input Num:1~256',1,p2);
% p = str2num(p3{1}); p = p/255;
%% 图片分割
image = im2bw(I, p); %小于阈值的为黑，大于阈值的为白
subplot(2, 2, 2); 
imshow(image); 
title('（b）图像前景与背景区分明显的分割结果')

2. 万有引力算法

2.1. 入口程序

function [Fbest,Lbest,BestChart,MeanChart]=GSA_image(objfun,N,max_it,ElitistCheck,min_flag,Rpower,...
    counts,x)
% 说明
% Main function for Gravitational Search Algorithm.
% V: 速度 Velocity.
% a: 加速度 Acceleration.
% M: 惯性质量 Mass. Ma=Mp=Mi=M;
% dim: 自变量维度 Dimension of the test function.
% N: 种群规模 Number of agents.
% X: 个体位置集，一个N*dim矩阵 Position of agents. dim-by-N matrix.
% R: 个体距离 Distance between agents in search space.
% [low-up]: 参数范围 Allowable range for search space.
% Rnorm: 范数，参考文献公式8 Norm in eq.8.
% Rpower: 参考文献公式7 Power of R in eq.7.

Rnorm=2; % 使用二阶范数 
% 获取目标函数参数界限、维数 get allowable range and dimension of the test function.
low=0.01;
up=255;
dim=1;%2
% 随机初始化种群 random initialization for agents.
X=initialization(dim,N,up,low); 
% 用于保存当前最优值和平均适应度值变化情况 create the best so far chart and average fitnesses chart.
BestChart=zeros(1,max_it);
MeanChart=zeros(1,max_it);
% 初始化个体解
V=zeros(N,dim);
for iteration=1:max_it
    % 检查是否越界 Checking allowable range. 
    X=space_bound(X,up,low);
    % 计算个体适应度函数值 Evaluation of agents.
    fitness=zeros(1,N);
    for agent=1:N
        fitness(1,agent)=objfun(X(agent,:),counts,x);
    end
    % 寻找当前迭代最优个体
    if min_flag==1
        [best,best_X]=min(fitness); % 最小化情况 minimization.
    else
        [best,best_X]=max(fitness); % 最大化情况 maximization.
    end
    if iteration==1
       Fbest=best;
       Lbest=X(best_X,:);
    end
    % 更新目前为止最优个体
    if min_flag==1
        if best<Fbest % 最小化情况 minimization.
            Fbest=best;
            Lbest=X(best_X,:);
        end
    else
        if best>Fbest % 最大化情况 maximization
            Fbest=best;
            Lbest=X(best_X,:);
        end
    end
    BestChart(iteration)=Fbest;
    MeanChart(iteration)=mean(fitness);
    % 计算惯性质量M（文献公式14—20） Calculation of M. eq.14-20
    M=massCalculation(fitness,min_flag);
    % 计算引力常亮（文献公式13） Calculation of Gravitational constant. eq.13.
    G=Gconstant(iteration,max_it);
    % 计算加速度 Calculation of accelaration in gravitational field. eq.7-10,21.
    a=Gfield(M,X,G,Rnorm,Rpower,ElitistCheck,iteration,max_it);
    % 个体移动 Agent movement. eq.11-12
    [X,V]=move(X,a,V);
    X
end

2.2 初始化种群程序

% This function initializes the position of the agents in the search space, randomly.
function X=initialization(dim,N,up,down)
if size(up,2)==1
    X=rand(N,dim).*(up-down)+down;
end
if size(up,2)>1
    for i=1:dim
        high=up(i);
        low=down(i);
        X(:,i)=rand(N,1).*(high-low)+low;
    end
end

2.3 检查是否越界

 %This function checks the search space boundaries for agents.
function X=space_bound(X,up,low)
[N,dim]=size(X);
for i=1:N
    % 对越界值进行重新初始化 Agents that go out of the search space, are reinitialized randomly.
    Tp=X(i,:)>up;
    Tm=X(i,:)<low;
    X(i,:)=(X(i,:).*(~(Tp+Tm)))+((rand(1,dim).*(up-low)+low).*logical((Tp+Tm)));
    % 将越界值重置为边界值 Agents that go out of the search space, are returned to the boundaries.
    % Tp=X(i,:)>up;
    % Tm=X(i,:)<low;
    % X(i,:)=(X(i,:).*(~(Tp+Tm)))+up.*Tp+low.*Tm;
end

2.4 对某一阈值适应度的计算程序

function f=objfun_image(cv,counts,x)
% cv为长度为2的横向量，即SVM中参数c和v的值

T=cv(1);

%% 选择图片，并二值化
% countsx=counts.*x;
sumI=sum(counts);
baifen=counts/sumI;
i=floor(T);
w0=sum(baifen(1:i));
w1=1-w0;
u0=sum(counts(1:i).*x(1:i))/sum(counts(1:i));
u1=sum(counts(i+1:length(x)).*x(i+1:length(x)))/sum(counts(i+1:length(x)));
f=w0*w1*(u0-u1)*(u0-u1);

2.5 计算惯性质量

% This function calculates the mass of each agent. eq.14-20
function M =massCalculation(fit,min_flag)
% Here, make your own function of 'mass calculation'
Fmax=max(fit);
Fmin=min(fit);
[~,N]=size(fit);
if Fmax==Fmin
    M=ones(N,1);
else
    if min_flag==1 % for minimization
        best=Fmin;
        worst=Fmax; % eq.17-18.
    else % for maximization
        best=Fmax;
        worst=Fmin; % eq.19-20.
    end    
    M=(fit-worst)./(best-worst); % eq.15.
end
M=M./sum(M); % eq.16.

2.6 计算引力常亮

% This function calculates Gravitational constant. eq.13.
function G=Gconstant(iteration,max_it)
% here, make your own function of 'G'.
alfa=20;
G0=100;
G=G0*exp(-alfa*iteration/max_it); % eq.28.

2.7 计算加速度程序

% This function calculates the accelaration of each agent in gravitational field. eq.7-10,21.
function a=Gfield(M,X,G,Rnorm,Rpower,ElitistCheck,iteration,max_it)
[N,dim]=size(X);
% In the last iteration, only 2 percent of agents apply force to the others.
% 在最后一次迭代中，只有百分之二的个体对其它个体有引力？？？
final_per=2; 
% 计算总引力 total force calculation
if ElitistCheck==1
    kbest=final_per+(1-iteration/max_it)*(100-final_per); % 参考文献公式21中kbest的计算 kbest in eq.21.
    kbest=round(N*kbest/100);
else
    kbest=N; % eq.9.
end
[~,ds]=sort(M,'descend');
E=zeros(N,dim);
for i=1:N % 遍历种群
    E(i,:)=zeros(1,dim);
    for ii=1:kbest
        j=ds(ii);
        if j~=i
            R=norm(X(i,:)-X(j,:),Rnorm); % 欧氏距离 Euclidian distanse.
            for k=1:dim
                E(i,k)=E(i,k)+rand*(M(j))*((X(j,k)-X(i,k))/(R^Rpower+eps));
                % note that Mp(i)/Mi(i)=1
            end
        end
    end
end
% 加速度 acceleration
a=E.*G; % note that Mp(i)/Mi(i)=1

2.8 计算个体移动程序

% This function updates the velocity and position of agents.
function [X,V]=move(X,a,V)
% movement.
[N,dim]=size(X);
V=rand(N,dim).*V+a; % eq.11.
X=X+V; % eq.12.

缺点

将迭代次数增加到20的时候会出现0/0导致的崩溃。

[1] 戚娜, 马占文. 基于万有引力搜索算法图像分割的实现[J]. 太赫兹科学与电子信息学报, 2017, 15(3):475-479.
[2] 齐丽娜, 张博, 王战凯. 最大类间方差法在图像处理中的应用[J]. 无线电工程, 2006, 36(7):25-26.
[3] 范炜锋. 万有引力搜索算法的分析与改进[D]. 广东工业大学, 2014,8-10.
[4] http://blog.csdn.net/u013337691/article/details/52732631