[fly,map] = imread('butterfly.gif'); % load image into MATLAB
figure(1)
image(fly),colormap(map); % display image
axis off, axis equal

fly=double(fly); % convert to double precision
[m,n]=size(fly);
mn = mean(fly,2); % compute row mean
X = fly - repmat(mn,1,n); % subtract row mean to obtain X
Z=1/sqrt(n-1)*X'; % create matrix, Z
covZ=Z'*Z; % covariance matrix of Z

%% Singular value decomposition
[U,S,V] = svd(covZ);
variances=diag(S).*diag(S); % compute variances
%bar(variances(1:30)); % scree plot of variances

%% Extract first 20 principal components
PCs=20;
VV=V(:,1:PCs);
Y=VV'*X; % project data onto PCs
ratio=256/(2*PCs+1) % compression ratio
XX=VV*Y; % convert back to original basis
XX=XX+repmat(mn,1,n); % add the row means back on
figure(2)
image(XX),colormap(map),axis off; % display results