PCA 图像融合 matlab+ python

可能存在错误内容 但是程序能通

可以从AI earth 网站下载不同波段的遥感图像数据（遥感图很大）
图像格式为tif 上传不了 用于融合的两张图像

RGB图像与灰度图像的PCA图像融合

% RGB图像与灰度图像的PCA图像融合，RGB得到通过PCA降维后三维数据，和灰度图像得到通过PCA降维后单维数据，然后直接进行替换，通过PCA逆变换(特征向量求逆后重构原矩阵)得到融合后的数据，有采用直方图配准的没有直接替换，不过估计直方图配准可以解决数据不同维度之间的问题，进行补数据
clear
clc
close all
img1 = double(imread('bldr_tm_6.tif'))/255; %归一化
img1=imresize(img1,[512,512]);% 设定图像大小
img3 = img1;
img2 = double(rgb2gray(imread('bldr_sp_1.tif')))/255;%全色图像 单光谱的 应当是单通道 但是由于图像读取图像是RGB三色 所以要转灰度
img2=imresize(img2,[512,512]);% 设定图像大小 保证两个图像降维后数据长度一致
figure;
% subplot(2,2,1);
imshow(img1,[]);title('多光谱图像');
figure;
% subplot(2,2,2);
imshow(img2,[]);title('全色光谱图像');

% PCA 过程 ,将各通道数据转为单维数据
[s1,s2,s3] = size(img1);
img1 = reshape(img1,[s1*s2,s3]); % s1*s2为图像长*宽 s3图像通道个数 重组成具有图像长*宽 长度的三个单维数据
% 零均值化（标准化） 这里标准化之后 数据会丢失
% ming = mean(img1);
% for i = 1:s3
% img1(:,i) = img1(:,i)-ming(:,i);% mean(img1)计算 img1三维度的均值
% end

% 求算协方差矩阵 由于上述已经减去过均值了 所以直接 各个变量与其他维度对应位置变量相乘 累加后的值 除以 单维变量总个数（图像像素个数） 就是 协方差矩阵
% covimg = (img1'*img1)/(s1*s2);% 两个矩阵的乘法 就是 各个变量与其他维度对应位置变量相乘 累加后的值，所以除以图像像素个数(size(img1,1)*size(img1,2))
covimg1 = sqrt((img1'*img1)/(s1*s2));% 开不开根号没区别
% 获取协方差矩阵特征值和特征向量
[V,D]=eig(covimg1);%[V,D] = eig(A) 返回特征值的对角矩阵 D 和特征向量矩阵 V，其列是对应的右特征向量，使得 A*V = V*D。
[d,ind] = sort(diag(D));%diag(D) 从 D 的对角线上提取特征值，然后按升序对得到的向量进行排序。sort 的第二个输出返回索引的置换向量。
Ds = D(ind,ind);%mathwork里面eig()函数使用有说明
Vs = V(:,ind);
% 获取各数据在所有成分上投影结果
for i =1 :s3
    img1vs(:,i) = img1*Vs(:,i);%将img1 所有 通道 的数据 投影到第 i 个主成分坐标轴上，其中i越大 对应主成分坐标的特征值越大，且 PCA逆变换为：img1vs(:,1)*Vs(:,i)’= img1 然后将所有的叠加即可还原
end

figure;
plot(img1vs(:,1),'.');
hold on 
plot(img1vs(:,2),'.');
hold on 
plot(img1vs(:,3),'.');
title('多光谱主成分降维数据')
legend(['第三主成分降维数据' ;'第二主成分降维数据' ;'第一主成分降维数据'])
% 获取 全色光谱的主成分数据 
% 由于全色光谱数据为单通道数据 一维度数据主成分是否为其本身（暂不确定，所以还是直接求一次）
[z1,z2] = size(img2);
img2 = reshape(img2,[z1*z2,1]); % s1*s2为图像长*宽 s3图像通道个数
covimg2 = sqrt((img2'*img2)/(z1*z2));
[V2,D2]=eig(covimg2);%[V,D] = eig(A) 返回特征值的对角矩阵 D 和矩阵 V，其列是对应的右特征向量，使得 A*V = V*D。
[d2,ind2] = sort(diag(D2));%diag(D) 从 D 的对角线上提取特征值，然后按升序对得到的向量进行排序。sort 的第二个输出返回索引的置换向量。
Ds2 = D2(ind2,ind2);%mathwork里面eig()函数使用有说明
Vs2 = V2(:,ind2);
img2vs = img2*Vs2;
figure;
plot(img2vs,'.');title('全色光谱降维数据')
% % 图像融合 就是将两个图像对应的主成分 数据进行替换 会有多种组合 此处选择第一主成分的数据进行替换
% 替换主成分数据 能实现图像融合 是因为主成分中包含了 色彩之间的差异信息 (PCA降维后保留数据之间的相对关系)
% 所以进行替换后能够提高色彩差异，并能够保留自身 的一些色彩差异
% PCA逆变换，将 替换后的主成分数据 进行 还原成图像 完成图像融合

n=3;
temp = img1vs(:,n); % n 表示要多光谱图像中要被替换的第n主成分 由于排序了 所以n越大对应的特征向量越大
img1vs(:,n) =0.3* img1vs(:,n) + 0.7*img2vs;% 0.3 和 0.7为权重系数

% subplot(2,2,3)
figure
imshow(reshape((img1vs)*(Vs'),[s1,s2,s3]),[]);title(['全色光谱降维数据 替换多光谱第 ' num2str(s3-n+1) ' 主成分数据']);% 由于n为升序 而第1主成分为最大，所以s3-n+1变一下
% subplot(2,2,4)
figure
imshow(reshape(temp*Vs2',[z1,z2]),[]);title([' 多光谱第 ' num2str(s3-n+1) ' 主成分数据 替换 全色光谱降维数据']);

% 被替换后的多光谱降维数据
figure;
plot(img1vs(:,1),'.');
hold on 
plot(img1vs(:,2),'.');
hold on 
plot(img1vs(:,3),'.');
title('替换后的多光谱主成分降维数据')
legend(['第三主成分降维数据' ;'第二主成分降维数据' ;'第一主成分降维数据'])

%%
% % 第一主成分提取得到色彩差异信息图像，
close all
img1vs1 = img1*Vs(:,end);
img1re1 = reshape(img1vs1,[s1,s2]);
% subplot(2,2,3);imshow(img1re1,[]);title('多光谱图像第一主成分投影结果图像');%是把多光谱中的色彩进行降维投影，第一主成分投影位置保留了最大色彩差异信息.

% %  第一主成分提取得到色彩差异信息图像，与 多光谱图像对应灰度图像的 差异
img1gray = rgb2gray(img3);
figure;
subplot(1,3,1);imshow(img1re1,[]);title('多光谱图像第一主成分投影结果图像');
subplot(1,3,2);imshow(img1gray,[]);title('多光谱图像的灰度图像')
sum(sum(abs(img1gray-img1re1)))%计算 多光谱图像的灰度图像 与 多光谱图像第一主成分投影结果图像之间的 差值 
subplot(1,3,3);
imshow(abs(img1gray-img1re1),[])% 显示 多光谱图像的灰度图像 与 多光谱图像第一主成分投影结果图像之间的 差值图像 （加了绝对值）
title('两者差异');

% %  最小特征值对应主成分提取得到色彩差异信息图像，及其与第一主成分提取得到色彩差异信息图像， 的差异
img1vs2 = img1*Vs(:,1);%最小特征值对应 主成分提取得到色彩差异信息图像 肯定是一个模糊图像
img1re2 = reshape(img1vs2,[s1,s2]);
figure;
sum(sum(abs(img1re2-img1re1)))
subplot(1,3,1);imshow(img1re1 ,[]);title('多光谱图像第一主成分投影结果图像')
subplot(1,3,2);imshow(img1re2 ,[]);title('多光谱图像最小特征值对应主成分投影结果图像')
subplot(1,3,3);imshow(abs(img1re1-img1re2) ,[]);title('两者差异图像')

%  所有主成分的投影结果
% img1vs = img1*Vs;
% img1re = reshape(img1vs,[s1,s2,s3]);
% figure;imshow(img1re,[]);title('多光谱图像所有主成分投影结果图像');%保留了所有色彩差异信息.一个投影维度（主成分）表示一个色彩差异信息分布，那这多个投影结果叠加好像没啥意义

% pca(img1);
% 
figure;
imshow(reshape(img1,[s1,s2,s3]),[])
figure;
img1vs3 = img1vs1*V(:,end)';
img1vs3 = reshape(img1vs3,[s1,s2,s3]);
imshow(img1vs3,[])

RGB1 和 RGB2图像融合

% RGB1 和 RGB2 分别重组成两个一维数据，2维数据放一起降维后
% 选取第一主成分投影数据作为图像融合结果 

clear
clc
close all
img2 = imread('bldr_tm_6.tif');
img1 = (imread('bldr_sp_1.tif'));
a = size(img1);
if ~isequal(size(img1), size(img2))
    img2 = imresize3(img2, [a(1) a(2) a(3)]);
end

% Concatenate the two images into a single matrix
img = [1.3*img2(:) img1(:)]; % 1.3 为参数 调整两个图像之间的权重
% Perform PCA on the image matrix

%这个不归一化 会失真严重 
[coeff, score, latent] = pca(double(img)./255);

% Create the fused image by taking the first principal component
fused_img = reshape(score(:, 1), size(img1));

% Display the fused image
figure;
subplot(1,3,1)
imshow(fused_img, []);title('融合图像')
subplot(1,3,2)
imshow(img1, []);title('img1')
subplot(1,3,3)
imshow(img2, []);title('img2')

% RGB1 和 RGB2 组成一个6维数据，6维放一起降维后
% 选取第一主成分投影数据作为图像融合结果 
clear
clc
close all
% 这个归一化和不归一化 没看出区别 不知道是不是因为是灰度图像 差异不大，而且在这归一化才行

img2 = double(imread('bldr_tm_6.tif'))./255;
img1 = double(imread('bldr_sp_1.tif'))./255;

% Resize the images to the same size if they have different dimensions
if ~isequal(size(img1), size(img2))
    img2 = imresize3(img2, size(img1));
end



% Concatenate the three images into a single matrix
[m,n,c]=size(img1);
img = cat(3, 0.8*img1, img2);% 0.8 为参数 调整两个图像之间的权重
img = reshape(img,[m*n,6]);

[coeff, score, latent] = pca(double(img));

figure;
for i=1:6
    subplot(2,3,i)
    imshow(reshape(score(:,i), [m,n]), []);title(['第' num2str(i) ' 主成分数据成像结果']);
end

figure;
imshow(img1, []);title('img1')
figure;
imshow(img2, []);title('img2')

补上python 版 （matlab 自带配准方法不好用 ）

import numpy as np
import cv2 as cv
import matplotlib.pyplot as plt

img1 = cv.cvtColor(cv.imread('2.tif'),cv.COLOR_BGR2RGB)  # 多光谱 RGB 
img2 = cv.cvtColor(cv.imread('1.tif'),cv.COLOR_BGR2RGB)  # 全色光谱 读出来也是RGB 
plt.figure(figsize=([10,10]))
plt.imshow(img2) # 显示全色 图像 
plt.show()

# sift 图像配准 
sift_detector = cv.SIFT_create()
# Find the keypoints and descriptors with SIFT
kp1, des1 = sift_detector.detectAndCompute(img1, None)
kp2, des2 = sift_detector.detectAndCompute(img2, None)
# BFMatcher with default params
bf = cv.BFMatcher()
matches = bf.knnMatch(des1, des2, k=2)
# Filter out poor matches
good_matches = []
for m,n in matches:
    if m.distance < 0.89*n.distance:
        good_matches.append(m)
matches = good_matches
points1 = np.zeros((len(matches), 2), dtype=np.float32)
points2 = np.zeros((len(matches), 2), dtype=np.float32)
for i, match in enumerate(matches):
    points1[i, :] = kp1[match.queryIdx].pt
    points2[i, :] = kp2[match.trainIdx].pt
# Find homography
H, mask = cv.findHomography(points1, points2, cv.RANSAC)
# Warp image 1 to align with image 2
img1Reg = cv.warpPerspective(img1, H, (img2.shape[1], img2.shape[0]))
# cv.imwrite('aligned_img1.jpg', img1Reg)#可以用于保存配准后的图像 
# The problem is that this matrix H is found via a compute-intensive optimization process
plt.figure(figsize=([10,10]))
plt.imshow(img1Reg)# 显示配准后的图像
plt.show()

from sklearn.decomposition import PCA  # 加载PCA算法包

img3 = np.reshape(img1Reg,( img1Reg.shape[0]*img1Reg.shape[1],3)) # 用np.reshape 重组 多光谱数据为RGB 转为 二维数据，第二维为3 
img4 = cv.cvtColor(img2,cv.COLOR_RGB2GRAY)  # 全色图转灰度
img4 = np.reshape(img4,(img4.shape[0]*img4.shape[1],1)) # 用np.reshape 重组 全色光谱数据为 二维数据，第二维为1

plt.figure(figsize=([10,10]))
plt.imshow(np.reshape(img5[:,0],(1390,1071)),'gray') #将多光谱数据为RGB 转为 二维数据，中的第一维数据转成图像
# plt.imshow(np.reshape(img4[:,0],(1390,1071)),'gray')
plt.show()

pca1 = PCA(n_components=3)  # 加载PCA算法，设置降维后主成分数目为3
reduced_x1 = pca1.fit_transform(img3)  # 对多光谱进行降维
pca2 = PCA(n_components=1)  # 加载PCA算法，设置降维后主成分数目为1
reduced_x2 = pca2.fit_transform(img4)  # 对全色进行降维

reduced_x1[:,0] = reduced_x2[:,0] + 1*reduced_x1[:,0] # 用全色降维数据 加权替换 多光谱降维数据  其中多光谱降维数据的0表示第一主成分数据 

reduced_x = np.dot(reduced_x1, pca1.components_) + pca1.mean_  # 还原数据，得到融合后的数据
reduced= np.reshape(reduced_x,(img2.shape[0],img2.shape[1],img2.shape[2]))# 重组融合后的数据成图像

plt.figure(figsize=([10,10]))
plt.imshow(np. uint8(reduced),'gray')# 显示重组融合后的图像
# plt.imshow(np.reshape(img4[:,0],(1390,1071)),'gray')
plt.show()

plt.figure(figsize=([10,10]))
plt.subplot(1,2,1)
plt.imshow(img1Reg,'gray')
plt.subplot(1,2,2)
plt.imshow(img2,'gray')
plt.show()