%=====================================================
% KLT demo for image processing course, 
% TELIN, Ghent University
%
% This demo is a modification of the "kltdemo"
% from DIP Spring 99 Class Demo
% Shreekanth Mandayam, Rowan University
%
% Aleksandra Pizurica, TELIN, Ghent University, 2003.
%=====================================================


clear; close all;
%loading images and displaying 
load test1.mat;
load test2.mat;
load test3.mat;
load test4.mat;
load test5.mat;
load test6.mat;


figure('Units','normalized','position',[0.01 0.5 0.9 0.25]);
subplot(1,5,1),imagesc(X1); colormap(gray);
title('Input image 1');
subplot(1,5,2),imagesc(X2); colormap(gray);
title('Input image 2');
subplot(1,5,3),imagesc(X3); colormap(gray);
title('Input image 3');
subplot(1,5,4),imagesc(X4); colormap(gray);
title('Input image 4');
subplot(1,5,5),imagesc(X5); colormap(gray);
title('Input image 5');


%Show al the components in one image
CompImg(:,1:50)=X1;
CompImg(:,50)=100;
CompImg(:,51:(51+49))=X2;
CompImg(:,100)=100;
CompImg(:,101:(101+49))=X3;
CompImg(:,150)=100;
CompImg(:,151:(151+49))=X4;
CompImg(:,200)=100;
CompImg(:,201:(201+49))=X5;
figure('Units','normalized','position',[0.01 0.2 0.9 0.25]);
imagesc(CompImg);
colormap(gray);
title('Input components');


[m n]=size(X1);
N=m*n;

%forming mean vector
meanx(1,1)=mean(mean(X1));
meanx(2,1)=mean(mean(X2));
meanx(3,1)=mean(mean(X3));
meanx(4,1)=mean(mean(X4));
meanx(5,1)=mean(mean(X5));
%meanx(6,1)=mean(mean(X6));

%forming covariance matrix;
p1=reshape(X1',1,N);
r1=reshape(X2',1,N);
d1=reshape(X3',1,N);
b1=reshape(X4',1,N);
t1=reshape(X5',1,N);
g1=reshape(X6',1,N);
X(1,:)=p1;
X(2,:)=r1;
X(3,:)=d1;
X(4,:)=b1;
X(5,:)=t1;
%X(6,:)=g1;

C=0;
for i=1:N;
 C=C+X(:,i)*(X(:,i))';
end;
C=C/N;
 
%computing eigenvalues and eigenvectors
[V,D]=eig(C);
D

%Computing Transform
for i=1:N;
 Y(:,i)=(V')*(X(:,i)-meanx);
end;

%Show the 3 Principal Transformed Images
im1=reshape(Y(5,:),m,n);
im2=reshape(Y(4,:),m,n);
im3=reshape(Y(3,:),m,n);
%im4=reshape(Y(2,:),m,n);

%Sorting and arranging
A(:,1)=V(:,5);
A(:,2)=V(:,4);
A(:,3)=V(:,3);
%A(:,4)=V(:,2);
%A=V;
A=A';

%Stored Ys correspond to these principal images
Ys(1,:)=Y(5,:);
Ys(2,:)=Y(4,:);
Ys(3,:)=Y(3,:);
%Ys(4,:)=Y(2,:);

%Ys=Y;

%Reconstruction
for i=1:N;
 Xr(:,i)=(A)'*Ys(:,i)+meanx;
end;




figure('Units','normalized','position',[0.01 0.2 0.6 0.25]);
subplot(1,3,1),imagesc(im1); colormap(gray);
title('Principal component 1');
subplot(1,3,2),imagesc(im2); colormap(gray);
title('Principal component 2');
subplot(1,3,3),imagesc(X3); colormap(gray);
title('Principal component 3');

