基于粒子群算法的图像配准程序
基于粒子群优化(PSO)算法的图像配准MATLAB实现。图像配准的目标是找到最优的空间变换参数，使两幅图像最佳对齐。
核心原理
图像配准的PSO方法将变换参数（如平移、旋转、缩放）作为粒子位置，通过优化相似性度量函数来找到最佳配准参数。
% 变换模型:
% [x'] [s*cosθ -s*sinθ tx] [x]
% [y'] = [s*sinθ s*cosθ ty] [y]
% [1 ] [0 0 1 ] [1]
完整PSO图像配准实现
1. 主程序框架
function pso_image_registration()
% 基于PSO的图像配准主程序
clear; clc; close all;
%% 1. 读取和预处理图像
fprintf('正在加载和预处理图像...\n');
[fixed_img, moving_img] = load_and_preprocess_images();
%% 2. PSO参数设置
pso_params = struct();
pso_params.n_particles = 20; % 粒子数量
pso_params.max_iter = 50; % 最大迭代次数
pso_params.w = 0.729; % 惯性权重
pso_params.c1 = 1.494; % 个体学习因子
pso_params.c2 = 1.494; % 社会学习因子
% 参数搜索范围 [tx, ty, θ, scale]
% tx, ty: 平移参数 (像素)
% θ: 旋转角度 (度)
% scale: 缩放因子
search_space = [-50, -50, -45, 0.5; % 下限
50, 50, 45, 2.0]; % 上限
%% 3. 执行PSO优化
fprintf('开始PSO优化...\n');
[best_params, best_fitness, convergence] = ...
pso_image_optimizer(fixed_img, moving_img, pso_params, search_space);
%% 4. 结果显示
fprintf('\n=== 配准结果 ===\n');
fprintf('最优变换参数:\n');
fprintf('平移 tx = %.2f 像素\n', best_params(1));
fprintf('平移 ty = %.2f 像素\n', best_params(2));
fprintf('旋转 θ = %.2f 度\n', best_params(3));
fprintf('缩放 scale = %.2f\n', best_params(4));
fprintf('相似性度量: %.6f\n', best_fitness);
%% 5. 可视化结果
visualize_registration_results(fixed_img, moving_img, best_params, convergence);
end
2. 图像加载和预处理
function [fixed_img, moving_img] = load_and_preprocess_images()
% 图像加载和预处理
% 方法1: 从文件加载（请修改为您的图像路径）
% fixed_img = imread('fixed_image.jpg');
% moving_img = imread('moving_image.jpg');
% 方法2: 生成测试图像
[fixed_img, moving_img] = create_test_images();
% 转换为灰度图像（如果是彩色图像）
if size(fixed_img, 3) == 3
fixed_img = rgb2gray(fixed_img);
end
if size(moving_img, 3) == 3
moving_img = rgb2gray(moving_img);
end
% 转换为double类型并归一化
fixed_img = im2double(fixed_img);
moving_img = im2double(moving_img);
% 显示原始图像
figure('Name', '原始图像', 'Position', [100, 100, 1200, 400]);
subplot(1,3,1);
imshow(fixed_img);
title('参考图像 (Fixed)');
subplot(1,3,2);
imshow(moving_img);
title('待配准图像 (Moving)');
end
function [fixed_img, moving_img] = create_test_images()
% 创建测试图像对
fixed_img = zeros(256, 256);
% 在参考图像中创建一些特征
fixed_img(80:120, 80:120) = 1; % 正方形
fixed_img(150:180, 150:200) = 0.7; % 矩形
fixed_img = imgaussfilt(fixed_img, 2); % 高斯模糊
% 对待配准图像应用变换（模拟配准问题）
tx_true = 15; % 真实平移x
ty_true = -10; % 真实平移y
theta_true = 10; % 真实旋转角度（度）
scale_true = 1.2; % 真实缩放因子
tform = affine2d([scale_true*cosd(theta_true) -scale_true*sind(theta_true) 0;
scale_true*sind(theta_true) scale_true*cosd(theta_true) 0;
tx_true ty_true 1]);
moving_img = imwarp(fixed_img, tform, 'OutputView', imref2d(size(fixed_img)));
% 添加噪声模拟真实情况
moving_img = imnoise(moving_img, 'gaussian', 0, 0.01);
end
3. PSO图像优化器
function [global_best, global_best_fitness, convergence] = ...
pso_image_optimizer(fixed_img, moving_img, pso_params, search_space)
n_particles = pso_params.n_particles;
dim = 4; % [tx, ty, θ, scale]
% 初始化粒子群
particles_pos = zeros(n_particles, dim);
particles_vel = zeros(n_particles, dim);
personal_best_pos = zeros(n_particles, dim);
personal_best_fitness = -inf(n_particles, 1); % 负无穷，因为我们要最大化相似性
% 在搜索空间内随机初始化位置和速度
for i = 1:n_particles
for j = 1:dim
range = search_space(2,j) - search_space(1,j);
particles_pos(i,j) = search_space(1,j) + rand() * range;
particles_vel(i,j) = -0.1 * range + 0.2 * range * rand();
end
end
% 全局最优初始化
global_best = particles_pos(1,:);
global_best_fitness = -inf;
convergence = zeros(pso_params.max_iter, 1);
%% PSO主循环
for iter = 1:pso_params.max_iter
% 动态惯性权重（线性递减）
w = pso_params.w * (1 - 0.5 * iter/pso_params.max_iter);
for i = 1:n_particles
% 计算当前粒子的适应度（相似性度量）
current_fitness = calculate_similarity(fixed_img, moving_img, particles_pos(i,:));
% 更新个体最优（最大化相似性）
if current_fitness > personal_best_fitness(i)
personal_best_fitness(i) = current_fitness;
personal_best_pos(i,:) = particles_pos(i,:);
end
% 更新全局最优
if current_fitness > global_best_fitness
global_best_fitness = current_fitness;
global_best = particles_pos(i,:);
end
end
% 更新所有粒子的速度和位置
for i = 1:n_particles
for j = 1:dim
r1 = rand();
r2 = rand();
% 速度更新
cognitive = pso_params.c1 * r1 * (personal_best_pos(i,j) - particles_pos(i,j));
social = pso_params.c2 * r2 * (global_best(j) - particles_pos(i,j));
particles_vel(i,j) = w * particles_vel(i,j) + cognitive + social;
% 位置更新
particles_pos(i,j) = particles_pos(i,j) + particles_vel(i,j);
% 边界约束处理
if particles_pos(i,j) < search_space(1,j)
particles_pos(i,j) = search_space(1,j);
particles_vel(i,j) = -0.5 * particles_vel(i,j);
elseif particles_pos(i,j) > search_space(2,j)
particles_pos(i,j) = search_space(2,j);
particles_vel(i,j) = -0.5 * particles_vel(i,j);
end
end
end
convergence(iter) = global_best_fitness;
% 显示进度
if mod(iter, 5) == 0
fprintf('迭代 %d/%d, 当前最优相似度: %.6f\n', iter, pso_params.max_iter, global_best_fitness);
end
end
end
4. 相似性度量函数
function similarity = calculate_similarity(fixed_img, moving_img, params)
% 计算变换后图像与参考图像的相似性
% params = [tx, ty, theta, scale]
tx = params(1);
ty = params(2);
theta = params(3);
scale = params(4);
% 创建仿射变换矩阵
tform = affine2d([scale*cosd(theta) -scale*sind(theta) 0;
scale*sind(theta) scale*cosd(theta) 0;
tx ty 1]);
try
% 应用变换
transformed_img = imwarp(moving_img, tform, 'OutputView', imref2d(size(fixed_img)));
% 计算相似性度量（多种选择）
% 1. 归一化互相关 (NCC) - 对亮度变化鲁棒
fixed_norm = fixed_img - mean(fixed_img(:));
transformed_norm = transformed_img - mean(transformed_img(:));
numerator = sum(fixed_norm(:) .* transformed_norm(:));
denominator = sqrt(sum(fixed_norm(:).^2) * sum(transformed_norm(:).^2));
if denominator == 0
similarity = 0;
else
similarity = numerator / denominator;
end
% 2. 互信息 (可选，计算量较大但更鲁棒)
% similarity = mutual_information(fixed_img, transformed_img);
% 3. 均方误差 (MSE) - 需要转换为最小化问题
% mse = mean((fixed_img(:) - transformed_img(:)).^2);
% similarity = 1 / (1 + mse); % 将MSE转换为相似性度量
catch
% 如果变换失败，返回很低的相似性
similarity = -inf;
end
end
function mi = mutual_information(img1, img2)
% 计算两幅图像的互信息
% 将图像量化为256级
img1_quant = round(img1 * 255) + 1;
img2_quant = round(img2 * 255) + 1;
% 计算联合直方图
joint_hist = accumarray([img1_quant(:), img2_quant(:)], 1, [256, 256]);
joint_prob = joint_hist / sum(joint_hist(:));
% 计算边缘概率
prob1 = sum(joint_prob, 2);
prob2 = sum(joint_prob, 1);
% 计算互信息
mi = 0;
for i = 1:256
for j = 1:256
if joint_prob(i,j) > 0 && prob1(i) > 0 && prob2(j) > 0
mi = mi + joint_prob(i,j) * log2(joint_prob(i,j) / (prob1(i) * prob2(j)));
end
end
end
end
5. 结果可视化
function visualize_registration_results(fixed_img, moving_img, best_params, convergence)
% 可视化配准结果
%% 应用最优变换
tx = best_params(1);
ty = best_params(2);
theta = best_params(3);
scale = best_params(4);
tform = affine2d([scale*cosd(theta) -scale*sind(theta) 0;
scale*sind(theta) scale*cosd(theta) 0;
tx ty 1]);
registered_img = imwarp(moving_img, tform, 'OutputView', imref2d(size(fixed_img)));
%% 创建可视化图形
figure('Name', '图像配准结果', 'Position', [100, 100, 1400, 800]);
% 1. 收敛曲线
subplot(2,3,1);
plot(1:length(convergence), convergence, 'b-', 'LineWidth', 2);
xlabel('迭代次数');
ylabel('相似性度量');
title('PSO收敛曲线');
grid on;
% 2. 原始图像对比
subplot(2,3,2);
imshow(fixed_img);
title('参考图像');
subplot(2,3,3);
imshow(moving_img);
title('待配准图像');
% 3. 配准后的图像
subplot(2,3,4);
imshow(registered_img);
title('配准后的图像');
% 4. 差异图像（配准前）
subplot(2,3,5);
diff_before = abs(fixed_img - moving_img);
imshow(diff_before, []);
title('配准前差异');
colorbar;
% 5. 差异图像（配准后）
subplot(2,3,6);
diff_after = abs(fixed_img - registered_img);
imshow(diff_after, []);
title('配准后差异');
colorbar;
%% 重叠显示
figure('Name', '图像重叠显示', 'Position', [200, 200, 1000, 400]);
% RGB合成显示
overlay = zeros([size(fixed_img), 3]);
overlay(:,:,1) = fixed_img; % 红色通道 - 参考图像
overlay(:,:,2) = registered_img; % 绿色通道 - 配准图像
overlay(:,:,3) = 0; % 蓝色通道
subplot(1,2,1);
imshow(overlay);
title('参考图像(红) vs 配准图像(绿)');
% 等高线叠加
subplot(1,2,2);
imshow(fixed_img);
hold on;
[C, h] = contour(registered_img, 5, 'r-', 'LineWidth', 1.5);
title('参考图像 + 配准图像等高线');
legend('配准图像轮廓');
%% 定量评估
fprintf('\n=== 定量评估 ===\n');
% 配准前的MSE
mse_before = mean((fixed_img(:) - moving_img(:)).^2);
fprintf('配准前MSE: %.6f\n', mse_before);
% 配准后的MSE
mse_after = mean((fixed_img(:) - registered_img(:)).^2);
fprintf('配准后MSE: %.6f\n', mse_after);
% 改进百分比
improvement = (mse_before - mse_after) / mse_before * 100;
fprintf('MSE改进: %.2f%%\n', improvement);
end
高级功能扩展
1. 多分辨率配准
function best_params = multi_resolution_pso(fixed_img, moving_img, pso_params, search_space, n_levels)
% 多分辨率PSO配准
best_params = [0, 0, 0, 1]; % 初始参数
for level = n_levels:-1:1
% 创建金字塔图像
scale_factor = 1 / (2^(level-1));
fixed_small = imresize(fixed_img, scale_factor);
moving_small = imresize(moving_img, scale_factor);
% 调整搜索范围
current_search_space = search_space;
current_search_space(1:2, 1:2) = current_search_space(1:2, 1:2) * scale_factor;
fprintf('正在处理金字塔层 %d (缩放因子: %.3f)\n', level, scale_factor);
% 在当前分辨率执行PSO
[best_params, ~] = pso_image_optimizer(fixed_small, moving_small, pso_params, current_search_space);
% 为下一层调整参数（平移参数需要缩放）
if level > 1
best_params(1:2) = best_params(1:2) / scale_factor;
end
end
end
2. 鲁棒相似性度量
function similarity = robust_similarity(fixed_img, transformed_img)
% 鲁棒相似性度量组合
% 1. 归一化互相关
ncc = calculate_ncc(fixed_img, transformed_img);
% 2. 梯度方向相似性
gradient_sim = calculate_gradient_similarity(fixed_img, transformed_img);
% 3. 结构相似性 (SSIM)
ssim_val = ssim(transformed_img, fixed_img);
% 加权组合
weights = [0.5, 0.3, 0.2];
similarity = weights(1) * ncc + weights(2) * gradient_sim + weights(3) * ssim_val;
end
function gs = calculate_gradient_similarity(img1, img2)
% 计算梯度方向相似性
[gx1, gy1] = gradient(img1);
[gx2, gy2] = gradient(img2);
mag1 = sqrt(gx1.^2 + gy1.^2);
mag2 = sqrt(gx2.^2 + gy2.^2);
% 梯度方向余弦相似性
dot_product = gx1.*gx2 + gy1.*gy2;
magnitudes = mag1 .* mag2;
valid_mask = magnitudes > 0;
if any(valid_mask(:))
gs = mean(dot_product(valid_mask) ./ magnitudes(valid_mask));
else
gs = 0;
end
end
参考代码 基于粒子群算法的图像配准程序 www.3dddown.com/cna/77861.html
建议和优化
性能优化：
- 图像金字塔：使用多分辨率策略加速收敛
- ROI选择：只在感兴趣区域计算相似性
- 并行计算：使用
parfor并行处理粒子评估
参数调优：
% 针对不同类型图像的推荐参数
medical_images_params = struct(...
'n_particles', 25, ...
'max_iter', 100, ...
'w', 0.7, ...
'c1', 1.5, ...
'c2', 1.5);
remote_sensing_params = struct(...
'n_particles', 30, ...
'max_iter', 80, ...
'w', 0.6, ...
'c1', 1.7, ...
'c2', 1.7);
浙公网安备 33010602011771号