% 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;

figure(3)
imshow(uint8(RGB));
title('PCA-RGB表示的轉化圖像');




              %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
              %                       下面是計算相關系數                           %
              %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
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));
fprintf('\n\n   R的相關系數為:%.4f\n   G的相關系數為:%.4f\n   B的相關系數為:%.4f',rou_R,rou_G,rou_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
fprintf('\n\n   R的清晰度為:%.4f\n   G的清晰度為:%.4f\n   B的清晰度為:%.4f',...
          g_R/(M-1)/(N-1),g_G/(M-1)/(N-1),g_B/(M-1)/(N-1));