一、擷取代碼方式
二、圖像融合與PCA算法簡介
1 圖像融合概念與分類
圖像融合是将針對同一場景不同來源的圖像通過特定的算法生成一張新的圖像，融合後的圖像綜合
了不同源圖像間的互補資訊，提高了後期圖像可被分析和加工處理的能力。圖像融合主要分為像素級融合、特征級融合以及決策級融合。
像素級融合直接處理原始像素點，要求不同源生成的圖像比對度高，是以圖像保留的資訊量大，但是
運算量巨大：特征級融合需要先提取各源圖像的特征資訊，然後利用權重的特征資訊進行圖像重建，該融合方式容易丢失圖像的部分資訊：決策級融合是綜合運用像素級融合與特征級融合，通過對圖像的分析依據具體情況找到最優比對。圖像融合分類示意見圖。
融合後的圖像一般需要進行品質評估，主觀方法就是人眼判斷融合圖像的品質好壞。客觀方法包含與
源圖像進行相關名額的計算、比對，例如均值、标準差、平均梯度、資訊熵、空間頻率等.
2 PCA算法原理
主成分分析技術PC A(principal components analysis) ， 也叫做主分量分析， 基本數學原理是K-L正交分解。K-L變換是通過原始信号集合構造正交基，信号與基進行積分運算後，每個基對應的系數表示原始信号在該基下的投影值。在連續信号領域，可通過信号的内在特征構造核函數，不同的變換核函數
對應不同的基：在離散領域，變換核函數變為變換矩陣，當變換矩陣為協方差矩陣時，離散K-L變換即為PCA。
設原始信号在本征空間， PC A可以了解為一種信号降維算法，通過協方差矩陣構造出一組基，使得
賦範子空間基信号的權重之和與本征空間誤差最小，即投影誤差最小。PC A變換後的每個基的系數即為本征空間在子空間下的投影值，此時可以選取有限個基的不同組合向原始信号通近I]，對于非稀疏矩陣，理論上隻需少量包含資訊量大的基投影，即可在相當程度上近似于本征空間。設離散随機信号為X，X為n行p列矩陣，見式(1)。
現以列為投影方向構造子空間的基(以行為方向同理)，X可改寫為式(2)。
j=1，2…p：x為包含n個元素的列向量。
設變換矩陣為A，X變換後的矩陣為F，則F見式(3)。
當按列方向投影時， 需保證n多p， 此時A為pxp方陣：同理， 當按行方向投影時需保證n≤p， A為nxn
方陣，A就是要構造的基回。此時式(3)可寫為式(4)。
1.3 變換矩陣的計算
通過1.2節的原理分析，隻要構造出滿足3個條件的變換矩陣A，就可以對信号X進行主成分分析，得到主成分矩陣F，并利用F子列的任意組合得到重構信号X。在高等代數裡，按列投影時，變換矩陣A就是信号X的協方差矩陣的特征向量按列組成的矩陣。同理，按行投影時，變換矩陣A就是信号X的協方差矩陣的特征向量按行組成的矩陣[1]，證明略，隻給出
2 圖像的融合
對于灰階圖像采用PCA算法時， 所得主成分矩陣的列與圖像的列長度相等，重構圖像會出現栅欄效
應， 原因是PC A算法是一種降維算法， 将二維圖像降維成一維向量了，是以在對圖像采用該算法時，一般采用多通道圖像， 最簡單的多通道圖像就是RGB彩色圖像。
圖像融合基本原理：通過PC A正變換分别将2個待融合的三通道彩色圖像，降維成3個單通道圖像，即3個主成分，用另一幅圖像的某個主成分替換掉該圖像的某個主成分， 然後進行PC A逆變換， 得到的重構圖像即為融合圖像。
具體計算過程如下所述。
1)設可見光彩色圖像為Xi，紅外光僞彩色圖像為X， 兩圖像的分辨率都為MXN， 通道數都為3， 分辨
率不同時可通過尺度縮放保持一緻，通道數不同時可通過低通道數圖像直接取代高通道數圖像的主成分。
2)設X在3個通道下的圖像分别為X1，Xi2，X 13， 即大小都是M×N的矩陣， 将其變為包含MXN個元素的列向量，記為X1，X12，X’13。
3)将X1，X12，X13組合成X。X=Xxi此時X為行數為MXN， 列數為3的矩陣， 依次帶入式(7一11)，得到轉換矩陣A。
4)通過式(4)得到主成分矩陣F，此時的F是行數為MXN， 列數為3的矩陣， 提取第1列， 将其轉化為MXN的矩陣， 則恰好對應X的第1主成分Fu，Fi就是圖像Xi的第1主成分圖像，第2主成分圖像、第3主成分圖像分别為F12，F13。同理，可以求出圖像X的3個主成分圖像為F21，F22，F23。
5)将Fi，Fz(ij=1，2，3)互換，并按式(6)作逆變換，即可求出融合圖像Y，不同的互換組合可以求出不同的Y， Xu， X 2； (i=1.2.3) 不僅可以是RGB通道下的分解， 也可以是HIS空間、YUV空間下的分解[14]，如此可以得到更多的Y。
三、部分源代碼
% function Pca = PCA(TM);
clear
g_R=0; %r清晰度描述
g_G=0; %g清晰度描述
g_B=0; %b清晰度描述
h_R=0; %熵的描述
h_G=0;
h_B=0;
fenzi_R=0;
fenzi_G=0;
fenzi_B=0;
fenmu_up_R=0;
fenmu_up_G=0;
fenmu_up_B=0;
fenmu_low_R=0;
fenmu_low_G=0;
fenmu_low_B=0;
init_up_R=[];
init_up_G=[];
init_up_B=[];
init_low_R=[];
init_low_G=[];
init_low_B=[];
up=imread('high.jpg'); %讀圖像
low=imread('low.jpg');
figure(1)
imshow(up); %讀RGB數值
title('PCA-RGB表示的高分辨率圖像');
figure(2)
imshow(low);
title('PCA-RGB表示的低分辨率圖像');
[up_R]=double(up(:,:,1));
[up_G]=double(up(:,:,2));
[up_B]=double(up(:,:,3));
[low_R]=double(low(:,:,1));
[low_G]=double(low(:,:,2));
[low_B]=double(low(:,:,3));
[M,N,color]=size(up);
up_Mx = 0;
low_Mx=0;
for i = 1 : M
for j = 1 : N
up_S = [up_R(i,j),up_G(i,j),up_B(i,j)]'; % 生成由R，G， B組成的三維列向量
up_Mx = up_Mx + up_S;
low_S = [low_R(i,j),low_G(i,j),low_B(i,j)]';
low_Mx = low_Mx + low_S;
end
end
up_Mx = up_Mx / (M*N); % 計算三維列向量的平均值
low_Mx = low_Mx / (M*N);
up_Cx = 0;
low_Cx=0;
for i = 1 : M
for j = 1 : N
up_S = [up_R(i,j),up_G(i,j),up_B(i,j)]';
up_Cx = up_Cx + up_S*up_S';
low_S = [low_R(i,j),low_G(i,j),low_B(i,j)]';
low_Cx = low_Cx + low_S*low_S';
end
end
up_Cx = up_Cx / (M * N)- up_Mx*up_Mx'; % 計算協方差矩陳
low_Cx = low_Cx / (M * N)- low_Mx*low_Mx';
[up_A,up_latent] = eigs(up_Cx); % 協方差矩陳的特征向量組成的矩陳----PCA變換的系數矩陳,特征值
[low_A,low_latent] = eigs(low_Cx);
for i = 1 : M
for j = 1 : N
up_X = [up_R(i,j),up_G(i,j),up_G(i,j)]'; % 生成由R，G， B組成的三維列
up_Y = up_A'*up_X; % 每個象素點進行PCA變換正變換
up_Y = up_Y';
up_R(i,j) = up_Y(1); % 高分辨率圖檔的第1主分量
up_G(i,j) = up_Y(2); % 高分辨率圖檔的第2主分量
up_B(i,j) = up_Y(3); % 高分辨率圖檔的第3主分量
low_X = [low_R(i,j),low_G(i,j),low_G(i,j)]';
low_Y = low_A'*low_X;
low_Y = low_Y';
low_R(i,j) = low_Y(1); % 低分辨率圖檔的第1主分量
low_G(i,j) = low_Y(2); % 低分辨率圖檔的第2主分量
low_B(i,j) = low_Y(3); % 低分辨率圖檔的第3主分量
end
end
for i = 1 : M
for j = 1 : N
up_Y = [up_R(i,j),up_G(i,j),up_B(i,j)]'; % 生成由R，G， B組成的三維列向量
up_X = up_A*up_Y; % 每個象素點進行PCA變換反變換
up_X = up_X';
up_r(i,j) = up_X(1);
up_g(i,j) = up_X(2);
up_b(i,j) = up_X(3);
low_Y = [up_R(i,j),low_G(i,j),low_B(i,j)]';
low_X = low_A*low_Y;
low_X = low_X';
low_r(i,j) = low_X(1);
low_g(i,j) = low_X(2);
low_b(i,j) = low_X(3);
end
end
%RGB(:,:,1)=up_r;
%RGB(:,:,2)=up_g;
%RGB(:,:,3)=up_b;
RGB(:,:,1)=low_r;
RGB(:,:,2)=low_g;
RGB(:,:,3)=low_b;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 下面是計算相關系數 %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
init_up_R=ones(M,N)*mean(up_R(:));
init_up_G=ones(M,N)*mean(up_G(:));
init_up_B=ones(M,N)*mean(up_B(:));
init_low_R=ones(M,N)*mean(low_R(:));
init_low_G=ones(M,N)*mean(low_G(:));
init_low_B=ones(M,N)*mean(low_B(:));
for i=1:M
for j=1:N
fenzi_R=fenzi_R+(up_R(i,j)-init_up_R(i,j))*(low_R(i,j)-init_low_R(i,j));
fenmu_up_R=fenmu_up_R+(up_R(i,j)-init_up_R(i,j))^2;
fenmu_low_R=fenmu_low_R+(low_R(i,j)-init_low_R(i,j))^2;
fenzi_G=fenzi_G+(up_R(i,j)-init_up_G(i,j))*(low_R(i,j)-init_low_G(i,j));
fenmu_up_G=fenmu_up_G+(up_R(i,j)-init_up_G(i,j))^2;
fenmu_low_G=fenmu_low_G+(low_R(i,j)-init_low_G(i,j))^2;
fenzi_B=fenzi_B+(up_R(i,j)-init_up_B(i,j))*(low_R(i,j)-init_low_B(i,j));
fenmu_up_B=fenmu_up_B+(up_R(i,j)-init_up_B(i,j))^2;
fenmu_low_B=fenmu_low_B+(low_R(i,j)-init_low_B(i,j))^2;
end
end
rou_R=fenzi_R/(sqrt(fenmu_up_R*fenmu_low_R));
rou_G=fenzi_G/(sqrt(fenmu_up_G*fenmu_low_G));
rou_B=fenzi_B/(sqrt(fenmu_up_B*fenmu_low_B));
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 下面是計算清晰度G %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for ii=1:M-1
for jj=1:N-1
g_R=g_R+sqrt((((low_r(ii+1,jj)-low_r(ii,jj))^2+(low_r(ii,jj+1)-low_r(ii,jj))^2))/2);
g_G=g_G+sqrt((((low_g(ii+1,jj)-low_g(ii,jj))^2+(low_g(ii,jj+1)-low_g(ii,jj))^2))/2);
g_B=g_B+sqrt((((low_b(ii+1,jj)-low_b(ii,jj))^2+(low_b(ii,jj+1)-low_b(ii,jj))^2))/2);
end
end 四、運作結果
五、matlab版本及參考文獻
1 matlab版本
2014a
2 參考文獻
[1] 蔡利梅.MATLAB圖像處理——理論、算法與執行個體分析[M].清華大學出版社，2020.
[2]楊丹,趙海濱,龍哲.MATLAB圖像處理執行個體詳解[M].清華大學出版社，2013.
[3]周品.MATLAB圖像處理與圖形使用者界面設計[M].清華大學出版社，2013.