matlab有向图分割算法,基于万有引力搜索算法图像分割的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=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,:)

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,:)

% 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