1. 算法框架设计
%% 核心流程
1. 图像预处理 → 2. PSO参数初始化 → 3. 适应度函数计算 → 4. 粒子群迭代 → 5. 阈值输出
2. PSO参数设置
% 粒子群参数
n_particles = 30; % 粒子数量
max_iter = 100; % 最大迭代次数
w = 0.7; % 惯性权重
c1 = 1.5; c2 = 1.5; % 学习因子
dim = 2; % 搜索维度(二阈值时为2)
lb = [0, 0]; % 下界
ub = [255, 255]; % 上界
3. 多适应度函数实现
3.1 类间方差(Otsu准则)
function fitness = otsu_fitness(thresholds, hist)T1 = thresholds(1); T2 = thresholds(2);w0 = sum(hist(1:T1)); w1 = sum(hist(T1+1:T2)); w2 = sum(hist(T2+1:end));mu0 = sum((1:T1)' .* hist(1:T1)) / w0;mu1 = sum((T1+1:T2)' .* hist(T1+1:T2)) / w1;mu2 = sum((T2+1:end)' .* hist(T2+1:end)) / w2;fitness = w0*w1*(mu0-mu1)^2 + w1*w2*(mu1-mu2)^2;
end
3.2 最大熵准则
function fitness = entropy_fitness(thresholds, hist)T1 = thresholds(1); T2 = thresholds(2);H0 = -sum((hist(1:T1)./sum(hist(1:T1))).*log2(hist(1:T1)./sum(hist(1:T1))));H1 = -sum((hist(T1+1:T2)./sum(hist(T1+1:T2))).*log2(hist(T1+1:T2)./sum(hist(T1+1:T2))));H2 = -sum((hist(T2+1:end)./sum(hist(T2+1:end))).*log2(hist(T2+1:end)./sum(hist(T2+1:end))));fitness = H0 + H1 + H2;
end
3.3 区域均匀性
function fitness = region_uniformity(img, thresholds)T1 = thresholds(1); T2 = thresholds(2);mask1 = img < T1; mask2 = img >= T1 & img < T2; mask3 = img >= T2;uni1 = std2(img(mask1)); uni2 = std2(img(mask2)); uni3 = std2(img(mask3));fitness = 1 / (uni1 + uni2 + uni3 + eps);
end
3.4 梯度信息融合
function fitness = gradient_fitness(img, thresholds)[Gx, Gy] = imgradientxy(img);grad_mag = sqrt(Gx.^2 + Gy.^2);T1 = thresholds(1); T2 = thresholds(2);mask1 = grad_mag < T1; mask2 = grad_mag >= T1 & grad_mag < T2; mask3 = grad_mag >= T2;fitness = sum(mask1(:)) + 0.5*sum(mask2(:)) + 2*sum(mask3(:));
end
4. PSO主循环实现
%% 初始化粒子群
particles = lb + (ub-lb) .* rand(n_particles, dim);
velocities = 0.1*(ub-lb) .* (2*rand(n_particles, dim) - 1);
pbest = particles; % 个体最优
gbest = particles(1,:); % 全局最优%% 适应度计算
fitness = arrayfun(@(i) otsu_fitness(particles(i,:), imhist(img)), 1:n_particles);%% 迭代优化
for iter = 1:max_iter% 更新速度r1 = rand(n_particles, dim); r2 = rand(n_particles, dim);velocities = w*velocities + c1*r1.*(pbest - particles) + c2*r2.*(gbest - particles);velocities = min(max(velocities, -abs(ub-lb)), abs(ub-lb)); % 速度限制% 更新位置particles = particles + velocities;particles = min(max(particles, lb), ub); % 边界处理% 计算新适应度new_fitness = arrayfun(@(i) otsu_fitness(particles(i,:), imhist(img)), 1:n_particles);% 更新个体最优update_idx = new_fitness < fitness;pbest(update_idx,:) = particles(update_idx,:);fitness(update_idx) = new_fitness(update_idx);% 更新全局最优[min_fit, min_idx] = min(fitness);if min_fit < otsu_fitness(gbest, imhist(img))gbest = particles(min_idx,:);end% 混沌扰动(防止早熟)if mod(iter,10) == 0particles = tent_map(particles);end
end
5. 多阈值扩展方法
5.1 多维PSO优化(三阈值示例)
dim = 3; % 三阈值
lb = [0, 0, 0]; ub = [255, 255, 255];
% 适应度函数改为多类间方差计算
5.2 协作学习策略
% 将高维问题分解为多个子问题
sub_swarm1 = particles(:,1:2); % 前两个阈值
sub_swarm2 = particles(:,3:end); % 第三个阈值
% 各子群独立优化后合并结果
参考代码 将基本粒子群用于阈值灰度图像分割,同时给出多种适应度函数 www.youwenfan.com/contentcnk/66080.html
6. 优化策略
-
混沌初始化:使用Logistic映射生成初始粒子群,提升全局搜索能力
function x = logistic_map(n, r=4)x = zeros(n,1);x(1) = rand();for i=2:nx(i) = r*x(i-1)*(1-x(i-1));end end -
动态参数调整:根据迭代次数自适应调整惯性权重
w = 0.9 - 0.5*(iter/max_iter); % 线性递减 -
GPU加速:利用CUDA并行计算适应度
gpu_img = gpuArray(img); fitness = arrayfun(@(i) otsu_fitness(particles(i,:), imhist(gpu_img)), 1:n_particles);