图像处理 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()

文章来源