K-L变换

附带的读取图像的函数

Brovey变换

function [Out]=Brovey(msi,pan)
[row1,col1,bands1]=size(msi);
[row2,col2]=size(pan);
rowscale=row2/row1;
colscale=col2/col1;
Out=zeros(row2,col2,3);
msi=double(msi);
pan=double(pan);
%对多波段影像进行重采样
ReampleMSI=zeros(row2,col2,bands1);
for i=1:bands1
    ReampleMSI(:,:,i)=BilineInter(msi(:,:,i),row2,col2);
end
%Brovey变换
for k=1:bands1-1
    for i=1:row2
        for j=1:col2
            Out(i,j,k)= pan(i,j)*ReampleMSI(i,j,k)/(ReampleMSI(i,j,1)+ReampleMSI(i,j,2)+ReampleMSI(i,j,3));
        end
    end
end

end


%%
function M=BilineInter(datar,p,q)
%datar为低分辨率，[p,q]为要变换的大小
[m,n,k]=size(datar);
A=zeros(m+2,n+2,k);
%M为上采用的结果矩阵
M=zeros(p,q,k);
for t=1:k
A(2:m+1,2:n+1,t)=datar(:,:,t);%得到各个波段的值
end
%下面就是根据双线性内插法进行插值
for i=1:p
    for j=1:q
        I=ceil(i*m/p);
        J=ceil(j*n/q);
        u=I-i*m/p;
        v=J-j*n/q;
        a=[A(I,J,:),A(I,J+1,:);A(I+1,J,:),A(I+1,J+1,:)];
        w=[(1-u)*(1-v),(1-u)*v;(1-v)*u,u*v];
        M(i,j,:)=sum(sum(a.*w));
    end
end 
end

K-L变换

function Out=KLFusion(msidata,pan)
[~,~,bands]=size(msidata);
[row,col]=size(pan);
msidata1=zeros(row,col,bands);%根据msidata的尺寸声明msidata1的尺寸
for i=1:bands
    msidata1(:,:,i)=BilineInter(msidata(:,:,i),row,col);
end
%进行PCA变换之前先进性双线性内插求得上采用之后的图像
Out=PCA(msidata1,pan);%对多光谱数据做主成分变换
end

function Out=PCA(msidata,pan)
[row,col,bands]=size(msidata);
%先声明一个方差协方差矩阵
Q=zeros(bands,bands);%先声明一个方差协方差矩阵
bandmeans=zeros(1,bands);
msidata=double(msidata);
for i=1:bands
    bandmeans(1,i)=MyMean(msidata(:,:,i));%顺便求得每个波段的平均值以便于后面进行协方差的解算
end

for m=1:bands
    for n=1:bands 
        %外面两层循环就代表协方差矩阵的尺寸
        if m==n
            %m=n就直接计算方差就可以了
            Q(m,n)=Variance(msidata(:,:,m),bandmeans(1,m));
        else
            Q(m,n)=Cov(msidata(:,:,m),msidata(:,:,n),bandmeans(1,m),bandmeans(1,n));
        end
    end
end
%求得方差协方差矩阵之后进行特征矩阵和特诊向量的求解
%求解特征值和特征矩阵
[V,D]=eig(Q);%其中V是特征列向量拼接成为一个矩阵，D中diag为特征值
%将特征值进行排序，然后再进行正变换，将位于第一位置的多光谱但波段影像用全色影像代替
eigenvalues=diag(D);
%将原始的特征向量进行对应的调整,这一步直接在sort函数中与特征值排序同步进行
%将特征值排序之后特征向量必须跟着一起排序
[Transform,MSI]=Sort(eigenvalues,V,msidata);%得到对应的变换矩阵用于逆变换
%先做主成分变换
% MSI=MSI*Trasnform;%row*col*4*4*4===大致认为是row*col*4
% --------------TODO,实际上这种想法是错误的，我的MSI中的每一个波段进行行列变换才可以进行相乘
ReshapeMSI=zeros(row*col,bands);%首先给ReshapeMSI预分配内存为了得到更好的运行速度

for i=1:bands
    ReshapeMSI(:,i) = reshape(MSI(:,:,i),[row*col,1]);%将每个波段转变为列向量
end
%做主成分变换
PCAPositiveTransform=ReshapeMSI*Transform;%(row2*col2,4)*(4,4)
%对全色影像按照 ReshapeMSI(:,1)进行直方图匹配之后再进行替换
firstComponent=reshape(ReshapeMSI(:,1),[row,col]);%将波段变为正常的二维矩阵
%求得映射之后的图像
panMatch=HistogramMatch(pan,firstComponent);%其中pan为待匹配图像，firstComponent为匹配的参照图像
%将该图像代替第一主成分进行假彩色显示
panMatch=reshape(panMatch,[row*col,1]);%先变换为列向量
%再进行插入
PCAPositiveTransform(:,1)=panMatch;
%PCA反变换------这里面不做替换就可以不进行变换
PCANegativeTransform=PCAPositiveTransform*inv(Transform);
%其实PCA最后的图像融合效果就是一个假彩色融合的效果
Out=zeros(row,col,3);
Out(:,:,1)=reshape(PCANegativeTransform(:,1),[row,col]);%第一主成分（列向量转化为图像矩阵）
Out(:,:,2)=reshape(PCANegativeTransform(:,2),[row,col]);%第二主成分（列向量转化为图像矩阵）
Out(:,:,3)=reshape(PCANegativeTransform(:,3),[row,col]);%第三主成分（列向量转化为图像矩阵）
for i=1:3
    Out(:,:,i)=LineExtension(Out(:,:,i));
end
end

%%
function cov=Cov(Band1,Band2,mean1,mean2)
%计算两个波段的协方差提取成为一个函数,传入均值的目的是防止均值的重复计算降低效率
[row,col]=size(Band1);
sum=0;
for i=1:row
    for j=1:col
        sum=sum+(Band1(i,j)-mean1)*(Band2(i,j)-mean2);
    end
end
cov=sum/row/col;
end

%%
function variance=Variance(Band,mean)
%其实虽然方差只是协方差的一个特殊的例子，但是为了提高效率还是单独提取出来比较好
[row,col]=size(Band);
sum=0;
for i=1:row
    for j=1:col
        sum=sum+(Band(i,j)-mean)^2;%MALTAB还是非常方便的，所有函数都给我们重载好了
    end
end
variance=sum/row/col;
end

%%
function result=MyMean(A)
[row,col] = size(A);
result=0;
for i=1:row
    for j=1:col
        result=result+A(i,j); %这里需要注意的问题最后的数值是否会超过限制？？导致结果不准确----但好像又不会
    end
end
result=result/row/col;
end

%%
function [V,MSI]=Sort(eigenvalues,V,MSI)
eigensize=size(eigenvalues);
% temp1=eigenvalues;
%进行排序的时候需要知道每个位置是由什么位置变换过来
for i=1:eigensize(1)-1
    for j=2:eigensize(1)
        k=j-1;
        temp=eigenvalues(j,1);%保存当前的值
        temp1=V(:,j);
        temp2=MSI(:,:,j);
         if temp>eigenvalues(k,1)
            while k>=1&&temp>eigenvalues(k,1)
                eigenvalues(k+1,1)=eigenvalues(k,1);%往后移动
                V(:,k+1)=V(:,k);%将列向量交换
                MSI(:,:,k+1)=MSI(:,:,k);
                k=k-1;
            end
        end
        %找到插入的位置为j+1
         eigenvalues(k+1,1)=temp;
         V(:,k+1)=temp1;
         MSI(:,:,k+1)=temp2;
    end
end
% sortdata=eigenvalues;
% %寻找第一主成分是哪一个波段
% for i=1:eigensize(1)
%     if temp1(i,1)==sortdata(eigensize(1),1)
%         index=i;
%         break;
%     end
% end
end


%%
%双线性内插
function M=BilineInter(datar,p,q)
%datar为低分辨率，[p,q]为要变换的大小
[m,n,k]=size(datar);
A=zeros(m+2,n+2,k);
%M为上采用的结果矩阵
M=zeros(p,q,k);
for t=1:k
A(2:m+1,2:n+1,t)=datar(:,:,t);%得到各个波段的值
end
%下面就是根据双线性内插法进行插值
for i=1:p
    for j=1:q
        I=ceil(i*m/p);
        J=ceil(j*n/q);
        u=I-i*m/p;
        v=J-j*n/q;
        a=[A(I,J,:),A(I,J+1,:);A(I+1,J,:),A(I+1,J+1,:)];
        w=[(1-u)*(1-v),(1-u)*v;(1-v)*u,u*v];
        M(i,j,:)=sum(sum(a.*w));
    end
end 
end

%%
function src1=HistogramMatch(src1,src2) %两个相同尺寸的矩阵进行匹配
%灰度图像的直方图匹配
%第一步先将uint16的转化为uint8的图像便于在计算机处理0-255
%对src1和src2先在0-255上拉伸再进行直方图匹配
src1=LineExtension(src1);
src2=LineExtension(src2);
grayscale=256;
dist1=zeros(grayscale,1);
dist2=zeros(grayscale,1);
dist1(:,1)=AP(src1(:,:),grayscale);
dist2(:,1)=AP(src2(:,:),grayscale);%计算累计概率
%第二步进行灰度值重分配
%这里我假设src2是参考影像，目的是把src1指定到src2上
src1(:,:)=RuleDN(dist1(:,1),dist2(:,1),src1(:,:));
%第三步进行匹配图像的显示

end

%计算每个灰度级的累计概率
function cumulative=AP(src,grayscale)
[row,col]= size(src);
%就C++而言尽量避免uint与其他的数进行操作，所以我选择直接转化为double再进行操作
A=double(src);
dist=zeros(grayscale,1);%由于我暂时不知道如何读取图片的位图，所以暂时默认是256
for i=1:row
    for j=1:col
        index=round(A(i,j))+1;
        try 
            dist(index,1)=dist(index,1)+1;%这个是浙大数据结构中讲过的一个知识点
        catch 
            warning('Problem using function!');
        end
    end
end
%累计完频率之后进行概率计算，但是实际上考虑到效率和后续直方图的匹配根本用不到概率，所以我选择直接返回dist
%但是实际上上一句话也很矛盾,转变为double之后就是double，但是还是可以减少不必要的算数操作
%下一步计算每个灰度级出现的概率
dist=dist/(row*col);
%计算累计的概率,在MATLAB中很奇怪的一件事就是任何变量都不需要提前进行声明就可以使用
cumulative(1,1)=dist(1,1);
for i=2:grayscale
    cumulative(i,1)=cumulative(i-1,1)+dist(i,1);
end
end

function src1=RuleDN(cum1,cum2,src1) %cum为单波段的累计概率,src1为对应波段的概率
[row,col]=size(cum1);
src1_size=size(src1);%取得该单波段的尺寸
Map=zeros(row,col);%先进行一个for循环求得对应的映射矩阵
for i=1:row
    %对于1-256的DN值建立一一的映射关系
    mapDN=-1;
    minDis=100000;
    for j=1:row
        if abs(cum1(i,1)-cum2(j,1))<minDis
            minDis=abs(cum1(i,1)-cum2(j,1));
            mapDN=j-1;%真实的DN值为0-255，但是下标是1-256
        end
    end
    Map(i,1)=mapDN;%表示DN值为i的映射为mapDN值
end
%上面建立好对应的映射关系之后就进行DN值的总体映射
for i=1:src1_size(1)
    for j=1:src1_size(2)
        index=src1(i,j);
        %index为真实的DN值，在映射表中寻找对应的映射值进行映射
        src1(i,j)=Map(index+1,1);
    end
end
%上面就做好了直方图的匹配了
end

HSI融合变换

function HSIFusion(msi,pan)
%第一步，先进行msi的去极化线性拉伸
%只对前面三个波段进行线性拉伸
[row,col]=size(pan);
RGB=zeros(row,col,3);
for i=1:3
    temp=BilineInter(msi(:,:,i),row,col);
    RGB(:,:,i)=LineExtension(temp);
end
%拉伸到0-255之后就可以进行HSI变换
HSI=RGB2HSI(RGB);
%将I成分使用全色影像进行替换
%但是替换之前需要转换到0-255
pan=LineExtension(pan);
HSI(:,:,3)=pan;
%替换之后进行HSI反变换
Out=HSI2RGB(HSI);
imshow(uint8(Out));
end

%下面是用于第三步和第六步的HSI<=>RGB的变换函数
%%
function Dist=RGB2HSI(M)
[row,col,bands]=size(M);
Dist=zeros(size(M));
M=double(M);
for i=1:row
    for j=1:col
        R=M(i,j,1);
        G=M(i,j,2);
        B=M(i,j,3);
        Numerator=0.5*(2*R-G-B);
        Denominator=sqrt((R-G)^2+(R-B)*(G-B));
        if R== G && G ==B
            theta = 0;
        else
            theta=acos(Numerator/Denominator);
        end
        %由于是HSI变换，所以H是1，S是2，I是3
        Dist(i,j,3)=(R+G+B)/3;
        Dist(i,j,2)=1-3/(R+G+B)*min([R G B]);
        if G<B
            Dist(i,j,1)=360-theta;
        else
            Dist(i,j,1)=theta;
        end
    end
end
end

%%
function Image=HSI2RGB(HSI)
[lines,samples] = size(HSI(:,:,1));
Image = zeros(lines,samples,3);
for i = 1:lines
    for j = 1:samples                        
        if HSI(i,j,1)>=0 && HSI(i,j,1)<120
            Image(i,j,3) = HSI(i,j,3)*(1-HSI(i,j,2));
            Image(i,j,1) = HSI(i,j,3)*(1+HSI(i,j,2)*cos(pi*HSI(i,j,1)/180)/cos(pi*(60-HSI(i,j,1))/180));
            Image(i,j,2) = 3*HSI(i,j,3)- Image(i,j,1)-Image(i,j,3);            
        end
    
        if HSI(i,j,1)>=120 && HSI(i,j,1)<240
            Image(i,j,1) = HSI(i,j,3)*(1-HSI(i,j,2));
            Image(i,j,2) = HSI(i,j,3)*(1+HSI(i,j,2)*cos(pi*(HSI(i,j,1)-120)/180)/cos(pi*(180-HSI(i,j,1))/180));
            Image(i,j,3) = 3*HSI(i,j,3)- Image(i,j,1)- Image(i,j,2);          
        end
        
        if HSI(i,j,1)>=240 && HSI(i,j,1)<360
            Image(i,j,2) = HSI(i,j,3)*(1-HSI(i,j,2));
            Image(i,j,3) = HSI(i,j,3)*(1+HSI(i,j,2)*cos(pi*(HSI(i,j,1)-240)/180)/cos(pi*(360-HSI(i,j,1))/180));          
            Image(i,j,1) = 3*HSI(i,j,3)- Image(i,j,2)- Image(i,j,3);
        end
    end
end
end

%%
%下面这个函数用于使用双线性内插法进行上采样
function M=BilineInter(datar,p,q)
%datar为低分辨率，[p,q]为要变换的大小
[m,n,k]=size(datar);
A=zeros(m+2,n+2,k);
%M为上采用的结果矩阵
M=zeros(p,q,k);
for t=1:k
A(2:m+1,2:n+1,t)=datar(:,:,t);%得到各个波段的值
end
%下面就是根据双线性内插法进行插值
for i=1:p
    for j=1:q
        I=ceil(i*m/p);
        J=ceil(j*n/q);
        u=I-i*m/p;
        v=J-j*n/q;
        a=[A(I,J,:),A(I,J+1,:);A(I+1,J,:),A(I+1,J+1,:)];
        w=[(1-u)*(1-v),(1-u)*v;(1-v)*u,u*v];
        M(i,j,:)=sum(sum(a.*w));
    end
end 
end