clc
clear all
% egt the image and show it
a = getimg(3);
imshow(a);
size(a)
shg

class(a)

% it is uint8, so turn it into double
A = double(a);
class(A)

% compute the left covariance matrix , but without removing the mean, for simplicity
X = A*A';
size(X)

% get the right covariance matrix as well
Y = A'*A;
size(Y)

% compute the PCA for 
[ux,wx,vx] =svd(X,0);
size(wx)

%show the diagonal
spy(wx>0.00001)
shg

% same for the right hand covariance matrix
[uy,wy,vy] = svd(Y,0);
spy(wx>0.00001);shg
spy(wy>0.00001);shg

% compare the eigenvalues
WX=diag(wx);
WY=diag(wy);
WX(1:5)
WY(1:5)

% show the spectrum
plot(WX);shg

% one can see more using the log
plot(log(WX));shg
plot(log(WX)'.-');shg
plot(log(WX),'.-');shg

% show the cumulative eigenvalue spectrum, normalized to the trace
sum(diag(X))
sum(diag(wx))
CX = cumsum(WX)/sum(WX)
plot(CX)
shg
ux(1:5,1)
uy(1:5,1)
vy(1,1:5)
uy(1,1:5)

% now let us represent the data matrixc (image) itself
[u,w,v] = svd(A,0);
W = diag(w);

% compare the eigenvalues to WX
W(1:5)
W(1:5).*W(1:5)

% the eigenvectors also match
ux(1:5,1)
u(1:5,1)

% show the cumulative spectrum
C = cumsum(W)/sum(W);
plot(C)

% let us rebuild the image from the SVD decomposition, using all vectors
B = u*w*v';
imshow(uint8(B));shg

% the difference is small
dA = abs(A-B);
max(max(dA))
max(max(A))
plot(C)
shg

% this is a function to set the diagonal to zero past the n+1'th element
function ww = trunc(w,n)
    ww=w;
    for i=n+1:size(w,1)
        ww(i,i) = 0;
    end
end


% keep only a 100 modes
w100 = trunc(w,100);
spy(w>0.00001)
figure(2);spy(w100>0.00001)
shg
B100 = u*w100*v';
imshow(uint8(B100));shg

% try 50
B50 = u*trunc(w,50)*v';
imshow(uint8(B50));
shg

% now 5, and this is too few
B5 = u*trunc(w,5)*v';
imshow(uint8(B5));shg

% what is the compression ratio?
(243*320)/(5*243+320*5+5)
imshow(uint8(u*trunc(w,10)*v'));
imshow(uint8(u*trunc(w,20)*v'));
shg

% try 20 and 30 modes
B20 = u*trunc(w,20)*v';
B5 = u*trunc(w,5)*v';
B30 = u*trunc(w,30)*v';
imshow(uint8(B30));shg
figure(2);imshow(a);

% look at the difference image (residual)
d5 = abs(A-B30);
imshow(uint(d5));
imshow(uint8(d5));
shg

% we need to amplify the difference to see it
imshow(uint8(10*d5));shg
imshow(uint8(10*abs(B10-A));shg
imshow(uint8(10*abs(B100-A))); shg
imshow(uint8(100*abs(B100-A))); shg
imshow(uint8(10*abs(B10-A))); shg
imshow(uint8(10*abs(B20-A))); shg
imshow(uint8(10*abs(B30-A))); shg