基于压缩系数的粒子群优化和重力搜索算法用于多级图像阈值化matlab源码

2023-01-07 11:37:53

一、基于压缩系数的粒子群优化和重力搜索算法用于多级图像阈值化

1.1粒子群算法

粒子群算法（ Particle Swarm Optimization, PSO）最早是由Eberhart和Kennedy于1995年提出，它的基本概念源于对鸟群觅食行为的研究。设想这样一个场景：一群鸟在随机搜寻食物，在这个区域里只有一块食物，所有的鸟都不知道食物在哪里，但是它们知道当前的位置离食物还有多远。最简单有效的策略？寻找鸟群中离食物最近的个体来进行搜素。PSO算法就从这种生物种群行为特性中得到启发并用于求解优化问题。

用一种粒子来模拟上述的鸟类个体，每个粒子可视为N维搜索空间中的一个搜索个体，粒子的当前位置即为对应优化问题的一个候选解，粒子的飞行过程即为该个体的搜索过程．粒子的飞行速度可根据粒子历史最优位置和种群历史最优位置进行动态调整．粒子仅具有两个属性：速度和位置，速度代表移动的快慢，位置代表移动的方向。每个粒子单独搜寻的最优解叫做个体极值，粒子群中最优的个体极值作为当前全局最优解。不断迭代，更新速度和位置。最终得到满足终止条件的最优解。

算法流程如下：

1、初始化

首先，我们设置最大迭代次数，目标函数的自变量个数，粒子的最大速度，位置信息为整个搜索空间，我们在速度区间和搜索空间上随机初始化速度和位置，设置粒子群规模为M，每个粒子随机初始化一个飞翔速度。

2、个体极值与全局最优解

定义适应度函数，个体极值为每个粒子找到的最优解，从这些最优解找到一个全局值，叫做本次全局最优解。与历史全局最优比较，进行更新。

3、更新速度和位置的公式

4、终止条件

（1）达到设定迭代次数；（2）代数之间的差值满足最小界限

以上就是最基本的一个标准PSO算法流程。和其它群智能算法一样，PSO算法在优化过程中，种群的多样性和算法的收敛速度之间始终存在着矛盾．对标准PSO算法的改进，无论是参数的选取、小生境技术的采用或是其他技术与PSO的融合，其目的都是希望在加强算法局部搜索能力的同时，保持种群的多样性，防止算法在快速收敛的同时出现早熟收敛。

1.2重力搜索算法

引力搜索算法(Gravitational Search Algorithm，GSA)是Esmat Rashedi等人在2009年提出的一种随机性启发式搜索算法，这种算法的灵感来自于牛顿的万有引力定律与运动定律：1.任意两个质点有通过连心线方向上的力相互吸引，该引力大小与它们质量的乘积成正比与它们距离的平方成反比。2.力使物体获得加速度。

在GSA中，质点被抽象成解空间的一个解，解之间存在一个相互吸引力，这个吸引力由解的质量与两个解之间的距离确定，质点的质量被抽象成解的评估函数值。在解空间中，每个解由其他解对其的吸引力获得加速度，质量更大(评估函数值更优)所提供的加速度更大，从而使解向更优解的方向移动。

看到这里，有些小伙伴会以为GSA是个与PSO差不多的算法，是的，GSA与PSO的外层框架是一致的，但是，它们最关键的粒子移动策略却是不同的：

1.在PSO中，粒子的移动只使用两个最佳位置pbest和gbest。但是在GSA中，智能体的方向是根据所有其他智能体或部分较优的智能体获得的总力来计算的。
2.PSO使用一种内存来更新速度（由于pbest和gbest）。然而，GSA是无内存的，只有智能体的当前位置在更新过程中起作用。
3.在PSO中，更新不考虑解之间的距离，而在GSA中，力与解之间的距离成反比。
4.PSO模拟了鸟类的社会行为，而GSA的灵感来自于物理现象。

GSA的主要过程如下：

1. 确定搜索空间。
2. 随机初始化个体种群。
3. 对种群进行适应度评价。
4. 更新引力常量G，更新种群中最好的个体 best与最差的个体worst，更新个体的质量。
5. 计算每个个体在不同方向上的总引力。
6. 计算个体的加速度，基于此更新个体的速度。
7. 根据速度更新个体在解空间中的位置。
8. 重复步骤3-7直到达到停止标准。

算法流程：

二、部分代码


                      
%              CPSOGSA for Multilevel Image Thresholding
% 

%Publication: Rather, S. A., & Bala, P. S. (2021). Constriction Coefficient Based Particle Swarm Optimization and Gravitational Search Algorithm for Multilevel Image Thresholding. Expert Systems, doi: 10.1111/exsy.12717, Wiley, SCIE (I.F = 2.587)%% 

                                                            
    
clear all
close all
clc

% Parameter initialization
%      I = imread('Aeroplane.tiff');
   I = imread('Cameraman.tiff');
 
  level = 5; %% Threshold = level-1 
% 
 N_PAR = level;                          %number of thresholds (number of levels-1) (dimensiones)
 dim = N_PAR;  
% 
 n = 15;                                  % Size of the swarm " no of objects " %%% Default (n = 15)
 Max_Iteration  = 300;                    % Maximum number of "iterations"      %%% Default (Max_Iteration  = 300)
% 
if size(I,3) == 1 %grayscale image
[n_countR, x_valueR] = imhist(I(:,:,1));
end
Nt = size(I,1) * size(I,2); 
 
% % Lmax indicated color segments 0 - 256

Lmax = 256;   %256 different maximum levels are considered in an image (i.e., 0 to 255)

for i = 1:Lmax
    if size(I,3) == 1  
        %grayscale image
        probR(i) = n_countR(i) / Nt;
    end
end
if size(I,3) == 1
    up = ones(n,dim) * Lmax;
    low = ones(n,dim);
end
 tic
 RunNo  = 1;   
    for k = [ 1 : RunNo ]  
       [CPSOGSA_bestit,CPSOGSA_bestF,CPSOGSA_Fit_bests]= CPSOGSA(I, Lmax, n,Max_Iteration,low,up,dim, level, probR);
       BestSolutions1(k) = CPSOGSA_bestF; 
 disp(['Run # ' , num2str(k),'::' 'Best estimates =',num2str(CPSOGSA_bestit)]);         % CPSOGSA
    end  
% /* Boxplot Analysis */
   figure
   boxplot([BestSolutions1'],{'CPSOGSA'});
   color = [([1 0 0])];
   h = findobj(gca,'Tag','Box'); 
   for j=1:length(h) 
   patch(get(h(j),'XData'),get(h(j),'YData'),color(j));
   end 
   title ('\fontsize{15}\bf Aeroplane (k=2)');
   % %  title ('\fontsize{15}\bf  Cameraman (k=2)');
   xlabel('\fontsize{15}\bf Algorithms');
   ylabel('\fontsize{15}\bf Best Fitness Values');
   box on
% % % 

% /* Graphical Analysis*/
figure
 plot(CPSOGSA_Fit_bests,'DisplayName','CPSOGSA','Color','b','LineStyle','-','LineWidth',3);
 disp( ['Time_CPSOGSA =', num2str(toc)]); 
 title ('\fontsize{15}\bf Aeroplane (k=2)'); % k=2,4,6,8,10
 % %  title ('\fontsize{15}\bf Cameraman (k=2)');
 xlabel('\fontsize{15}\bf Iterations');
 ylabel('\fontsize{15}\bf Fitness values');
 legend('\fontsize{12}\bf CPSOGSA');
 %
 %
 gBestR = sort(CPSOGSA_bestit);
 Iout = imageGRAY(I,gBestR);
 Iout2 = mat2gray(Iout); 
 
% % Show results on images  

figure
imshow(Iout)
    
figure
imshow(I)
    
% % Show results

intensity = gBestR(1:dim-1);  
STDR  = std(CPSOGSA_Fit_bests)              %Standard deviation of fitness values       
MSEV = MSE(I, Iout)                         %Mean Square Error
PSNRV = PSNR(I, Iout)                       %PSNR between original image I and the segmented image Iout
SSIMV = ssim (I, Iout)                      %SSIM Quality Measure
FSIMV = FeatureSIM (I, Iout)                %FSIM Quality Measure
Best_Fitness_Value= CPSOGSA_Fit_bests(k)    %Best fitness
     
% % Plot the threshold values over the histogram
figure 
plot(probR)
    hold on
vmax = max(probR);
for i = 1:dim-1
    line([intensity(i), intensity(i)],[0 vmax],[1 1],'Color','r','Marker','.','LineStyle','-')
    title ('\fontsize{15}\bf Aeroplane (k=2)');
% title ('\fontsize{15}\bf  Cameraman (k=2)');
    xlabel('\fontsize{15}\bf Gray level');
    ylabel('\fontsize{15}\bf Frequency');
    hold off
end