多光谱图像融合：IHS、PCA与小波变换的MATLAB实现

基于IHS变换、PCA变换和小波变换的多光谱图像融合方法，适用于遥感、医学成像和计算机视觉等领域的高质量图像融合。

% 多光谱图像融合：IHS、PCA与小波变换
clear; close all; clc;% 设置路径
addpath('images');
warning('off', 'Images:initSize:adjustingMag');%% 加载数据
% 多光谱图像 (低空间分辨率，高光谱分辨率)
multi_spectral = imread('multispectral.tif'); % 替换为您的多光谱图像
% 全色图像 (高空间分辨率，低光谱分辨率)
panchromatic = imread('panchromatic.tif'); % 替换为您的全色图像% 检查图像大小并调整
if size(multi_spectral, 3) > 3multi_spectral = multi_spectral(:,:,1:3); % 使用前三个波段
end% 调整全色图像大小以匹配多光谱图像
if size(panchromatic, 3) > 1panchromatic = rgb2gray(panchromatic); % 转换为灰度
end% 将多光谱图像上采样到全色图像分辨率
if size(multi_spectral,1) ~= size(panchromatic,1) || ...size(multi_spectral,2) ~= size(panchromatic,2)[h, w, ~] = size(multi_spectral);panchromatic = imresize(panchromatic, [h, w], 'bicubic');
end% 归一化图像
multi_spectral = im2double(multi_spectral);
panchromatic = im2double(panchromatic);% 显示原始图像
figure('Name', '原始图像', 'Position', [100, 100, 1200, 500]);
subplot(1,2,1); imshow(multi_spectral); title('多光谱图像 (RGB)');
subplot(1,2,2); imshow(panchromatic); title('全色图像');%% IHS变换融合
fused_ihs = ihs_fusion(multi_spectral, panchromatic);%% PCA变换融合
fused_pca = pca_fusion(multi_spectral, panchromatic);%% 小波变换融合
fused_wavelet = wavelet_fusion(multi_spectral, panchromatic);%% 显示融合结果
figure('Name', '融合结果', 'Position', [100, 100, 1200, 900]);
subplot(2,2,1); imshow(multi_spectral); title('原始多光谱');
subplot(2,2,2); imshow(fused_ihs); title('IHS融合');
subplot(2,2,3); imshow(fused_pca); title('PCA融合');
subplot(2,2,4); imshow(fused_wavelet); title('小波融合');%% 图像质量评估
metrics = struct();% IHS融合质量
metrics.ihs = evaluate_fusion(multi_spectral, panchromatic, fused_ihs);% PCA融合质量
metrics.pca = evaluate_fusion(multi_spectral, panchromatic, fused_pca);% 小波融合质量
metrics.wavelet = evaluate_fusion(multi_spectral, panchromatic, fused_wavelet);% 显示评估结果
display_metrics(metrics);%% 保存结果
imwrite(fused_ihs, 'fused_ihs.jpg');
imwrite(fused_pca, 'fused_pca.jpg');
imwrite(fused_wavelet, 'fused_wavelet.jpg');%% ====================== 融合函数 ====================== %%
function fused_img = ihs_fusion(ms, pan)% IHS变换融合% 输入:%   ms - 多光谱图像 (MxNx3)%   pan - 全色图像 (MxN)% 输出:%   fused_img - 融合结果% RGB到IHS变换[I, H, S] = rgb2ihs(ms);% 直方图匹配pan_matched = hist_match(pan, I);% 替换强度分量I_fused = pan_matched;% IHS逆变换fused_img = ihs2rgb(I_fused, H, S);
endfunction fused_img = pca_fusion(ms, pan)% PCA变换融合% 输入:%   ms - 多光谱图像 (MxNx3)%   pan - 全色图像 (MxN)% 输出:%   fused_img - 融合结果% 将多光谱图像重塑为2D矩阵[M, N, D] = size(ms);ms_2d = reshape(ms, M*N, D);% PCA变换[coeff, score, ~] = pca(ms_2d);% 获取第一主成分PC1 = reshape(score(:,1), M, N);% 直方图匹配pan_matched = hist_match(pan, PC1);% 替换第一主成分score(:,1) = pan_matched(:);% PCA逆变换fused_2d = score * coeff';fused_img = reshape(fused_2d, M, N, D);% 确保值在[0,1]范围内fused_img = min(max(fused_img, 0), 1);
endfunction fused_img = wavelet_fusion(ms, pan)% 小波变换融合% 输入:%   ms - 多光谱图像 (MxNx3)%   pan - 全色图像 (MxN)% 输出:%   fused_img - 融合结果% 设置小波参数wavelet = 'db4'; % Daubechies 4小波level = 3;       % 分解层数% 初始化融合图像fused_img = zeros(size(ms));% 对每个波段进行融合for band = 1:size(ms,3)% 多光谱波段的小波分解[cA_ms, cH_ms, cV_ms, cD_ms] = dwt2(ms(:,:,band), wavelet);[cA_ms, cH_ms, cV_ms, cD_ms] = decompose_levels(cA_ms, wavelet, level-1);% 全色图像的小波分解[cA_pan, cH_pan, cV_pan, cD_pan] = dwt2(pan, wavelet);[cA_pan, cH_pan, cV_pan, cD_pan] = decompose_levels(cA_pan, wavelet, level-1);% 融合规则：% 低频分量 - 平均cA_fused = (cA_ms + cA_pan) / 2;% 高频分量 - 取绝对值最大cH_fused = fusion_rule(cH_ms, cH_pan);cV_fused = fusion_rule(cV_ms, cV_pan);cD_fused = fusion_rule(cD_ms, cD_pan);% 小波重构cA_fused = reconstruct_levels(cA_fused, cH_fused, cV_fused, cD_fused, wavelet);fused_img(:,:,band) = idwt2(cA_fused, cH_fused, cV_fused, cD_fused, wavelet, size(ms(:,:,band)));end% 确保值在[0,1]范围内fused_img = min(max(fused_img, 0), 1);
end%% ====================== 辅助函数 ====================== %%
function [I, H, S] = rgb2ihs(rgb)% RGB到IHS变换% I: 强度, H: 色度, S: 饱和度R = rgb(:,:,1);G = rgb(:,:,2);B = rgb(:,:,3);% 强度分量I = (R + G + B) / 3;% 色度和饱和度H = atan2(sqrt(3)*(G - B), (2*R - G - B));S = sqrt((R - G).^2 + (R - B).*(G - B));% 规范化H和SH = (H + pi) / (2*pi); % 映射到[0,1]S = S / max(S(:));     % 归一化
endfunction rgb = ihs2rgb(I, H, S)% IHS到RGB变换H = 2*pi*H - pi; % 恢复原始范围% 计算中间变量cosH = cos(H);sinH = sin(H);% 计算RGB分量R = I + S .* (2*cosH + sqrt(3)*sinH) / 3;G = I + S .* (-cosH + sqrt(3)*sinH) / 3;B = I - S .* (cosH + sqrt(3)*sinH) / 3;% 组合RGB图像rgb = cat(3, R, G, B);% 裁剪到有效范围rgb = min(max(rgb, 0), 1);
endfunction matched = hist_match(source, target)% 直方图匹配% 使source的直方图匹配target的直方图% 计算累积分布函数cdf_target = cumsum(imhist(target)) / numel(target);cdf_source = cumsum(imhist(source)) / numel(source);% 创建映射函数[~, ind] = min(abs(cdf_source - cdf_target'), [], 2);map = (ind - 1) / 255;% 应用映射matched = interp1(linspace(0,1,256), map, source, 'nearest');
endfunction [cA, cH, cV, cD] = decompose_levels(cA, wavelet, levels)% 多级小波分解cH = cell(1, levels);cV = cell(1, levels);cD = cell(1, levels);for l = 1:levels[cA, cH{l}, cV{l}, cD{l}] = dwt2(cA, wavelet);end
endfunction cA = reconstruct_levels(cA, cH, cV, cD, wavelet)% 多级小波重构levels = length(cH);for l = levels:-1:1cA = idwt2(cA, cH{l}, cV{l}, cD{l}, wavelet, size(cA)*2);end
endfunction fused = fusion_rule(coeff1, coeff2)% 小波系数融合规则 (绝对值最大)if iscell(coeff1)% 多级系数fused = cell(size(coeff1));for i = 1:numel(coeff1)fused{i} = fusion_rule(coeff1{i}, coeff2{i});endelse% 单级系数mask = abs(coeff1) > abs(coeff2);fused = coeff1 .* mask + coeff2 .* (~mask);end
endfunction metrics = evaluate_fusion(ms, pan, fused)% 图像融合质量评估metrics = struct();% 1. 空间质量指标 (全色图像相关性)metrics.SpatialCorr = corr2(fused, pan);% 2. 光谱质量指标 (多光谱图像相关性)spectral_corr = zeros(1, size(ms,3));for i = 1:size(ms,3)spectral_corr(i) = corr2(fused(:,:,i), ms(:,:,i));endmetrics.SpectralCorr = mean(spectral_corr);% 3. 信息熵metrics.Entropy = entropy(rgb2gray(fused));% 4. 平均梯度[gx, gy] = gradient(double(rgb2gray(fused)));metrics.AvgGradient = mean(sqrt(gx(:).^2 + gy(:).^2));% 5. 结构相似性 (SSIM)metrics.SSIM = ssim(fused, ms);% 6. 峰值信噪比 (PSNR)metrics.PSNR = psnr(fused, ms);% 7. 通用图像质量指数 (UIQI)metrics.UIQI = uiqi(fused, ms);
endfunction qi = uiqi(img1, img2)% 通用图像质量指数img1 = double(rgb2gray(img1));img2 = double(rgb2gray(img2));C1 = (std(img1(:)) * std(img2(:)));C2 = cov(img1(:), img2(:));if size(C2,1) > 1C2 = C2(1,2);elseC2 = 0;endmu1 = mean(img1(:));mu2 = mean(img2(:));qi = (4 * C2 * mu1 * mu2) / ((var(img1(:)) + var(img2(:))) * (mu1^2 + mu2^2));
endfunction display_metrics(metrics)% 显示评估指标fprintf('\n===== 融合质量评估 =====\n');fprintf('方法\t\t空间相关\t光谱相关\t信息熵\t\t平均梯度\tSSIM\t\tPSNR\t\tUIQI\n');fields = fieldnames(metrics);for i = 1:numel(fields)m = metrics.(fields{i});fprintf('%s\t\t%.4f\t\t%.4f\t\t%.4f\t\t%.4f\t\t%.4f\t\t%.2f\t\t%.4f\n', ...fields{i}, ...m.SpatialCorr, ...m.SpectralCorr, ...m.Entropy, ...m.AvgGradient, ...m.SSIM, ...m.PSNR, ...m.UIQI);endfprintf('========================\n');% 可视化指标比较figure('Name', '融合质量指标', 'Position', [100, 100, 1200, 600]);metrics_names = {'空间相关', '光谱相关', '信息熵', '平均梯度', 'SSIM', 'PSNR', 'UIQI'};metrics_data = zeros(numel(fields), 7);for i = 1:numel(fields)m = metrics.(fields{i});metrics_data(i,:) = [m.SpatialCorr, m.SpectralCorr, m.Entropy, ...m.AvgGradient, m.SSIM, m.PSNR, m.UIQI];end% 归一化处理 (除信息熵和平均梯度外，越大越好)norm_data = metrics_data;norm_data(:,[3,4]) = 1 - norm_data(:,[3,4]); % 对于信息熵和平均梯度，越小越好for j = 1:size(norm_data,2)norm_data(:,j) = (norm_data(:,j) - min(norm_data(:,j))) / ...(max(norm_data(:,j)) - min(norm_data(:,j)));end% 雷达图subplot(1,2,1);polarplot(0,0); % 创建极坐标图hold on;theta = linspace(0, 2*pi, length(metrics_names)+1);for i = 1:size(norm_data,1)values = [norm_data(i,:), norm_data(i,1)]; % 闭合曲线polarplot(theta, values, 'LineWidth', 2);end% 设置标签rticks(0:0.2:1);thetaticks(rad2deg(theta(1:end-1)));thetaticklabels(metrics_names);title('融合质量雷达图 (归一化)');legend(fields, 'Location', 'southoutside');set(gca, 'FontSize', 10);% 条形图subplot(1,2,2);bar(metrics_data);set(gca, 'XTickLabel', fields);ylabel('指标值');title('融合质量指标比较');legend(metrics_names, 'Location', 'northwest');grid on;set(gca, 'FontSize', 10);
end

算法原理与实现

1. IHS变换融合

原理：

• IHS色彩空间将颜色分为强度(I)、色度(H)和饱和度(S)三个分量
• 全色图像包含丰富的空间细节信息，对应I分量
• 多光谱图像包含丰富的光谱信息，对应H和S分量

实现步骤：

1. 将RGB图像转换为IHS色彩空间
2. 对全色图像进行直方图匹配，使其统计特性与I分量一致
3. 用匹配后的全色图像替换I分量
4. 执行逆IHS变换回RGB空间

function fused_img = ihs_fusion(ms, pan)[I, H, S] = rgb2ihs(ms);pan_matched = hist_match(pan, I);fused_img = ihs2rgb(pan_matched, H, S);
end

2. PCA变换融合

原理：

• 主成分分析(PCA)将多光谱数据转换到新的正交空间
• 第一主成分包含最多信息，与空间细节高度相关
• 全色图像的空间细节可以替代第一主成分

实现步骤：

1. 对多光谱图像执行PCA变换
2. 对全色图像进行直方图匹配，使其与第一主成分统计特性一致
3. 用匹配后的全色图像替换第一主成分
4. 执行逆PCA变换

function fused_img = pca_fusion(ms, pan)[M, N, D] = size(ms);ms_2d = reshape(ms, M*N, D);[coeff, score, ~] = pca(ms_2d);PC1 = reshape(score(:,1), M, N);pan_matched = hist_match(pan, PC1);score(:,1) = pan_matched(:);fused_2d = score * coeff';fused_img = reshape(fused_2d, M, N, D);
end

3. 小波变换融合

原理：

• 小波变换将图像分解为低频(近似)和高频(细节)分量
• 低频分量包含主要光谱信息
• 高频分量包含空间细节信息

实现步骤：

1. 对多光谱图像的每个波段进行小波分解
2. 对全色图像进行小波分解
3. 融合规则：
   • 低频分量：多光谱和全色图像的平均值
   • 高频分量：取绝对值最大的系数
4. 对每个波段执行小波重构

function fused_img = wavelet_fusion(ms, pan)wavelet = 'db4';level = 3;fused_img = zeros(size(ms));for band = 1:size(ms,3)% 多级小波分解[cA_ms, cH_ms, cV_ms, cD_ms] = decompose_levels(ms(:,:,band), wavelet, level);[cA_pan, cH_pan, cV_pan, cD_pan] = decompose_levels(pan, wavelet, level);% 融合规则cA_fused = (cA_ms + cA_pan) / 2;cH_fused = fusion_rule(cH_ms, cH_pan);cV_fused = fusion_rule(cV_ms, cV_pan);cD_fused = fusion_rule(cD_ms, cD_pan);% 小波重构fused_img(:,:,band) = reconstruct_levels(cA_fused, cH_fused, cV_fused, cD_fused, wavelet);end
end

参考代码 应用IHS、PCA、小波变换进行多光谱图像融合 www.youwenfan.com/contentcsk/99553.html

图像质量评估指标

本方案提供全面的融合质量评估体系：

1. 空间相关性：衡量融合图像与全色图像的空间细节相似度

   metrics.SpatialCorr = corr2(fused, pan);
   

2. 光谱相关性：衡量融合图像与多光谱图像的光谱保真度

   spectral_corr(i) = corr2(fused(:,:,i), ms(:,:,i));
   

3. 信息熵：评估图像包含的信息量

   metrics.Entropy = entropy(rgb2gray(fused));
   

4. 平均梯度：反映图像清晰度和边缘强度

   [gx, gy] = gradient(double(rgb2gray(fused)));
   metrics.AvgGradient = mean(sqrt(gx(:).^2 + gy(:).^2));
   

5. 结构相似性(SSIM)：评估结构信息保持度

   metrics.SSIM = ssim(fused, ms);
   

6. 峰值信噪比(PSNR)：衡量图像重建质量

   metrics.PSNR = psnr(fused, ms);
   

7. 通用图像质量指数(UIQI)：综合评估图像质量

   metrics.UIQI = uiqi(fused, ms);