clear;
clc;
tic;
% b=imread('E:\pic\图片\图片\baboon.png');
% b=imread('E:\实验结果\膨胀加噪实验结果\模糊核聚类ls-svm预测\残差\0014nsctjlsvm.bmp');
% b=imread('E:\pic\加噪\luoding.bmp');
b=imread('E:\图片\1-10.bmp');
% b=rgb2gray(b);
a=b(:,:,1);
% imhist(a);
[mz,nz]=size(a);
% 计算平均邻域灰度的一维灰度直方图
a0=double(a);
a1=zeros(mz+2,nz+2);
a1(2:mz+1,2:nz+1)=a0;
a1(1,2:nz+1)=a1(2,2:nz+1);
a1(mz+2,2:nz+1)=a1(mz+1,2:nz+1);
a1(2:mz+1,nz+2)=a1(2:mz+1,nz+1);
a1(2:mz+1,1)=a1(2:mz+1,2);
a1([1 mz+2],[1 nz+2])=a1([2 mz+1],[2 nz+1]);
a3=uint8(colfilt(a1,[3 3],'sliding',@mean));
a2=a3(2:mz+1,2:nz+1);
%计算二维直方图
p=zeros(256,256);
for i=1:mz
for j=1:nz
n1=a(i,j);
n2=a2(i,j);
p(n1+1,n2+1)=p(n1+1,n2+1)+1;
end
end
pp=p./mz./nz;
calp=zeros(511,1);
h=zeros(511,1);
for i=1:256
for j=1:256
if pp(i,j)~=0
calp(i+j-1)=calp(i+j-1)+pp(i,j);
h(i+j-1)=h(i+j-1)+pp(i,j)*log(pp(i,j));
end
end
end
aa=find(calp~=0);
start=aa(1);
last=aa(end);
%确定基本粒子群算法的迭代次数和粒子个数
Iter = 50;
Agents = 20;
N = 4;
%初始化基本粒子群优化算法的参数
vm=5;
Wmax = 0.9; %最大惯性因子
Wmin = 0.1; %最小惯性因子
c1 = 2.0; %个体学习因子
c2 = 2.0; %群体学习因子
xmin = start; %粒子允许的最小位置
ymin = start;
xmax = last-1; %粒子允许的最大位置
ymax = last-1;
Vxmin = -vm; %粒子允许的最小速度
Vymin = -vm;
Vxmax = vm; %粒子允许的最大速度
Vymax = vm;
Xmin = [xmin,ymin];
Xmax = [xmax,ymax];
Vmin = [Vxmin,Vymin];
Vmax = [Vxmax,Vymax];
%R = 2*sqrt(Lmax^2+Hmax^2);
R = 10*sqrt(2);
%初始化第一个种群，Xpos为位置向量(起点坐标x，y)，Vpos为速度向量
%位置向量初始化,速度向量初始化
x=round(unifrnd(xmin,xmax,Agents,1));
vx=round(unifrnd(Vxmin,Vxmax,Agents,1));
% y=round(unifrnd(ymin,ymax,Agents,1));
for i=1:Agents
y(i)=round(unifrnd(x(i)+1,ymax,1,1));
end
vy=round(unifrnd(Vymin,Vymax,Agents,1));
Xpop=[x,y',vx,vy];
%生成小生境种群
pop = zeros(Agents/N,4,N);
for i = 1:N
pop(1:Agents/N,1:4,i) = Xpop((i-1)*Agents/N+1:i*Agents/N,1:4);
end
%初始化粒子适应度pFit、历史最佳适应度pBestFit和全局最佳适应度gBestFit
pFit = zeros(N,Agents/N);
pBestFit = pFit;
gBestFit = zeros(N,1);
ind = zeros(N,1);
PregBestFit = 0;
%初始化粒子历史最优位置pBest和全局最优位置gBest
pBest = zeros(Agents/N,2,N);
gBest = zeros(N,2);
%初始化粒子停滞次数nc和适应度变化值VarFit
% nc = zeros(N,Agents/N);
% VarFit =
% zeros(N,Agents/N);??????????????????????????????????????????????????????
% calp=calobjvalue(I);
%迭代更新
for it = 1:Iter
%计算适应度pFit并更新历史最优适应度pBestFit和全局最佳适应度gBestFit
for j = 1:N
for t = 1:Agents/N
Theta = pop(t,1:2,j);
%计算当前粒子的适应度pFit
pFit(j,t) = ThetaFit3(Theta,calp,h);
%更新个体历史最优适应度pBestFit
if pFit(j,t) > pBestFit(j,t)
pBestFit(j,t) = pFit(j,t);
pBest(t,1:2,j) = pop(t,1:2,j);
end
end
end
%更新当前时刻子群全局最优位置gBest
%和子群全局最佳适应度gBestFit
for j = 1:N
[gBestFit(j),ind(j)] = max(pBestFit(j,:));
gBest(j,1:2) = pBest(ind(j),1:2,j);
end
%更新前一时刻全局最佳适应度函数值PregBestFit和全局最优位置PregBest
[PregBestFit,ind1] = max(gBestFit(:));
PregBest = gBest(ind1,1:2);
%RCS小生境进化
s =0;
for i = 1:N-1
for j = i+1:N
D = sqrt((gBest(i,1)-gBest(j,1))^2+(gBest(i,2)-gBest(j,2))^2);
if D < R
if gBestFit(i) < gBestFit(j)
gBestFit(i) = 0;
else
gBestFit(j) = 0;
end
end
end
end
for i = 1:N
if gBestFit(i) == 0 %对置零的最优个体重新初始化
gBest(i,1) = round(unifrnd(xmin,xmax));
gBest(i,2) = round(unifrnd(gBest(i,1)+1,ymax));
pop(ind(i),1:2,i) = gBest(i,1:2);
Theta = pop(ind(i),1:2,i);
pFit(i,ind(i)) = ThetaFit3(Theta,calp,h);
pBestFit(i,ind(i)) = pFit(i,ind(i));
pBest(ind(i),1:2,i) = pop(ind(i),1:2,i);
[gBestFit(i),ind(i)] = max(pBestFit(i,:));
gBest(i,1:2) = pBest(ind(i),1:2,i);
s = 1;
end
end
while(s ~= 0)
s = 0;
for i = 1:N-1
for j = i+1:N
D = sqrt((gBest(i,1)-gBest(j,1))^2+(gBest(i,2)-gBest(j,2))^2);
if D < R
if gBestFit(i) < gBestFit(j)
gBestFit(i) = 0;
else gBestFit(j) = 0;
end
end
end
end
for i = 1:N
if gBestFit(i) == 0 %对置零的最优个体重新初始化
gBest(i,1) = round(unifrnd(xmin,xmax));
gBest(i,2) = round(unifrnd(gBest(i,1)+1,ymax));
pop(ind(i),1:2,i) = gBest(i,1:2);
Theta = pop(ind(i),1:2,i);
pFit(i,ind(i)) = ThetaFit3(Theta,calp,h);
pBestFit(i,ind(i)) = pFit(i,ind(i));
pBest(ind(i),1:2,i) = pop(ind(i),1:2,i);
[gBestFit(i),ind(i)] = max(pBestFit(i,:));
gBest(i,1:2) = pBest(ind(i),1:2,i);
s = 1;
end
end
end
%对最劣粒子重新初始化
if mod(it,10) == 0 %达到一定代数
[G,Ind] = min(gBestFit(:));
for i = 1:Agents/N
pop(i,1,Ind) = round(unifrnd(xmin,xmax));
pop(i,2,Ind) = round(unifrnd(pop(i,1,Ind)+1,ymax));
end
for j = 1:Agents/N
Theta = pop(j,1:2,Ind);
pFit(Ind,j) = ThetaFit3(Theta,calp,h);
if pFit(Ind,j) > pBestFit(Ind,j)
pBestFit(Ind,j) = pFit(Ind,j);
pBest(j,1:2,Ind) = pop(j,1:2,Ind);
end
end
[gBestFit(Ind),ind(Ind)] = max(pBestFit(Ind,:));
gBest(Ind,1:2) = pBest(ind(Ind),1:2,Ind);
end
%变尺度混沌变异
b = 0;
Pc = zeros(2,1);
it0 = 0;
m = 0;
n = 2;
gBest1 = zeros(1,2);
for i = 1:N
gBest1 = gBest(i,1:2);
while(b==0||b==0.25||b==0.5||b==0.75)
b = rand;
end
while(it0 < 10)
b = 4*b*(1-b);
for j = 1:2
Pc(j) = Xmin(j)+b*(Xmax(j)-Xmin(j));
m = 1-((it-1)/it)^n;
gBest1(j) = round((1-m)*gBest1(j)+m*Pc(j));
end
Theta = gBest1(1:2);
Fit = ThetaFit3(Theta,calp,h);
if Fit > gBestFit(i)
gBest(i,1:2) = gBest1;
break;
end
it0 = it0+1;
end
end
%位置向量Xpos和速度向量Vpos更新
for i = 1:N
for j = 1:Agents/N
for t = 3:4
%惯性权重更新并产生随机数r1、r2
W = Wmax - it*(Wmax - Wmin)/Iter;
r1 = rand;
r2 = rand;
%速度向量Vpos更新
pop(j,t,i) = W*pop(j,t,i) + c1*r1*(pBest(j,t-2,i) - pop(j,t-2,i)) + ...
c2*r2*(gBest(i,t-2) - pop(j,t-2,i));
%限制更新后粒子速度在[Vmin,Vmax]范围内
if pop(j,t,i) > Vmax(t-2)
pop(j,t,i) = Vmax(t-2);
end
if pop(j,t,i) < Vmin(t-2)
pop(j,t,i) = Vmin(t-2);
end
%粒子位置更新
pop(j,t-2,i) = round(pop(j,t-2,i) + pop(j,t,i));
if pop(j,t-2,i) > Xmax(t-2)
pop(j,t-2,i) = Xmax(t-2);
end
if pop(j,t-2,i) < Xmin(t-2)
pop(j,t-2,i) = Xmin(t-2);
