% Script: ShowLowRank
% Shows that the typical 8x8 pixel patch in an image is well approximated
% by a short sum of rank-1 matrices that are based on the 8x8 DCT matrix.

close all
% Read in an image, convert to black and white, and display...
A = imread('GoldenGate.jpg');
A = rgb2gray(A);
imshow(A)
% Click on the picture and pull out an 8x8 patch...
axis ij
[j,i] = ginput(1);
clc
Patch = A(i:i+7,j:j+7)
% Draw the cross-hair...
A(i,:) = 0; A(:,j) = 0;
imshow(A)
% Magnify and display the patch...
figure
set(gcf,'position',[900 600 200 200])
A1 = uint8(kron(double(Patch),ones(30,30)));
imshow(A1)
set(gcf,'color','y')

% Now examine low-rank approximations of the patch using the DCT basis.
% Set up C8 and its inverse..
C8 = cos((pi/16)*((0:7)'*(1:2:15)));
IC8 = C8'; IC8(:,1) = IC8(:,1)/8;IC8(:,2:8) = IC8(:,2:8)/4;

% The C8 transform (as in the first step in JPEG compression)...
B = C8*double(Patch)*C8';
disp('C8*Patch*C8'' = ')
disp(' ')
for k=1:8
disp(sprintf('%8.1f',B(k,:)))
end

% Note that
% Patch = IC8*B*IC8' = sum_over_all_ij B(i,j)*IC8(:,i)*IC8(:,j)'
%
% That is, Patch is a linear combination of the 64 rank-1 matrices
% IC8(:,i)*IC8(:,j)' where i and j each range from 1 to 8.
%
% Let's build up this sum of rank-1 matrices, pausing to display what the
% approximation looks like as an image. Since the B(i,j) tend to be
% large in the upper left corner, we will add in those more important
% rank-1 matrices first. In particular, at the kth step we will add in
% B(i,j)IC8(:,i)*IC8(:,j)' where i = idx(k) and j = jdx(k) with

idx = [1 2 1 3 2 1 4 3 2 1 5 4 3 2 1 6 5 4 3 2 1 7 6 5 4 3 2 1 8 7 6 5 4 3 2 1];
jdx = [1 1 2 1 2 3 1 2 3 4 1 2 3 4 5 1 2 3 4 5 6 1 2 3 4 5 6 7 1 2 3 4 5 6 7 8];
idx = [idx 8 7 6 5 4 3 2 8 7 6 5 4 3 8 7 6 5 4 8 7 6 5 8 7 6 8 7 8];
jdx = [jdx 2 3 4 5 6 7 8 3 4 5 6 7 8 4 5 6 7 8 5 6 7 8 6 7 8 7 8 8];
% Initialize the running sum...
P = zeros(8,8);
figure
set(gcf,'position',[1300 600 200 200])
for k=1:64
i = idx(k);
j = jdx(k);
P = P + B(i,j)*IC8(:,i)*IC8(:,j)';
err(k) = max(max(abs(double(Patch)-P)));
A1 = uint8(kron(P,ones(30,30)));
imshow(A1)
set(gcf,'color','c')
title(sprintf('k = %1d B(%1d,%1d) = %6.2f err = %5.2f',k,i,j,B(i,j),err(k)),'fontsize',12)
shg
pause(.1)
end
figure
plot(err)
title('Max pixel difference for rank-k approximation')