處理效果：

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

總結

以上是生活随笔為你收集整理的matlab有向图分割算法,基于万有引力搜索算法图像分割的MATLAB实现的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。