【阙值分割】基于matlab粒子群算法自适应多阈值图像分割【含Matlab源码 1459期】

原创

Matlab科研 2021-11-08 14:19:22 ©著作权

文章标签 matlab 算法图像处理 matlab图像处理权重 文章分类 数据结构与算法人工智能

©著作权归作者所有：来自51CTO博客作者Matlab科研的原创作品，请联系作者获取转载授权，否则将追究法律责任

一、粒子群算法自适应多阈值图像分割简介

理论知识参考：【基础教程】基于matlab图像处理图像分割【含Matlab源码 191期】

粒子群优化的多阈值图像自分割算法

二、部分源代码

clc;clear;close all;                

%% 输入图像；
Imag = imread('24063.jpg');%296059   
Imag=rgb2gray(Imag);
Image_OSTU=Imag;

%% 开始种群等基本定义
N = 500;                           % 初始种群个数
d = 2;                             % 阈值个数（参看上述的函数表达式）
ger = 300;                         % 最大迭代次数     
plimit = [1,256];
          
vlimit = [-2.5, 2.5;-2.5, 2.5];    % 设置速度限制
w = 0.8;                           % 惯性权重,个体历史成绩对现在的影响0.5~1之间
%还有自适应调整权重、随机权重等等
%（不同的权重设置很影响性能，按需要选取）

c1 = 0.5;                          % 自我学习因子
c2 = 0.5;                          % 群体学习因子 

x = zeros(N, 2);  
for i = 1:N                                             %对每一个个体
    x(i,1) = floor(plimit(1) + (plimit(2) - plimit(1)) * rand);%初始种群的位置
    x(i,2) = floor(x(i,1) + (plimit(2) - x(i,1)) * rand);%初始种群的位置                         
end 
                                    
v = rand(N, d);                   % 初始种群的速度,500行2列分别在两个维度上
xm = x;                           % 每个个体的历史最佳位置
ym = zeros(1, d);                 % 种群的历史最佳位置，两个维度，设置为0
fxm = zeros(N, 1);                % 每个个体的历史最佳适应度，设置为0
fym = -inf;                       % 种群历史最佳适应度,求最大值先设置成负无穷


iter=1;                             %初始的迭代次数因为用while设置为一
times = 1; 
record = zeros(ger, 1);             %记录器

%% 迭代更新开始
while iter <= ger

fx = f(x);               % 代入x中的二维数据，算出个体当前适应度,为500行1列的数据   
     for i = 1:N                    %对每一个个体做判断
        if fxm(i) < fx(i)           %如果每个个体的历史最佳适应度小于个体当前适应度
            fxm(i) = fx(i);         % 更新个体历史最佳适应度,第一轮就是把小于零的清除
            xm(i,:) = x(i,:);       % 更新个体历史最佳位置
        end 
     end
     
if fym < max(fxm)                  %种群历史最佳适应度小于个体里面最佳适应度的最大值
        [fym, nmax] = max(fxm);    % 更新群体历史最佳适应度,取出最大适应度的值和所在行数即位置
      
end

 v = v * w + c1 * rand *(xm - x) + c2 * rand *(repmat(ym, N, 1) - x); 
                                   % 速度更新公式,repmat函数把ym矩阵扩充成N行1列

 %%边界速度处理

        for j=1:N
        if  v(j,i)>vlimit(i,2)      %如果速度大于边界速度，则把速度拉回边界
            v(j,i)=vlimit(i,2);
        end
        if  v(j,i) < vlimit(i,1)     %如果速度小于边界速度，则把速度拉回边界
            v(j,i)=vlimit(i,1);
        end
        end
 end      
   
 x = floor(x + v);                          % 位置更新
 
 
 
 for j=1:N
        if  x(j,1)> plimit(2)
            x(j,1)=plimit(2);
        end
        if  x(j,1) < plimit(1)
            x(j,1)=plimit(1);
        end
        if  x(j,2)> plimit(2)
            x(j,2)=plimit(2);
        end
        if  x(j,2) < x(j,1)
            x(j,2)=x(j,1);
        end
 end
  

end

%% 作图
figure(3);
plot(record);                %画出最大值的变化过程
title('收敛过程')

threshold1 = ym(1);
threshold2 = ym(2);
[height,length]=size(Image_OSTU);
for i=1:length
    for j=1:height
        if Image_OSTU(j,i)>=threshold2
            Image_OSTU(j,i)=255;
        elseif Image_OSTU(j,i)<=threshold1
            Image_OSTU(j,i)=0;
        else  
            Image_OSTU(j,i)=(threshold1+threshold2)/2;
        end
    end
end
figure(4);
imshow(Image_OSTU);
xlabel(['最大类间差法阈值',num2str(ym)]);


%% 适应度函数
function fx = f(x)
Imag = imread('24063.jpg');%296059   
Imag=rgb2gray(Imag);

[height,length]=size(Imag);
totalNum=height*length;

pixelCount=zeros(1,256);%统计各个像素值的个数
for i=1:length
    for j=1:height
        number=Imag(j,i)+1;
        pixelCount(number)=pixelCount(number)+1;
    end
end


a=1:256;

fx = zeros(1, 500); 
for i=1:500
    m=x(i,1);
    n=x(i,2);
    w0=sum(pi(1:m));
    w1=sum(pi(m+1:n));
    w2=sum(pi(n+1:256));
    mean0=sum(pi(1:m).*a(1:m))/w0;
    mean1=sum(pi(m+1:n).*a(m+1:n))/w1;
    mean2=sum(pi(n+1:256).*a(n+1:256))/w2;
    fx(i)=w0*w1*(mean0-mean1)^2+w0*w2*(mean0-mean2)^2+w1*w2*(mean1-mean2)^2;
end
end