clear;
clc;
close all;

files = dir('E:\Berkeley\images\train\*.jpg'); % download the Berkeley database from http://www.eecs.berkeley.edu/Research/Projects/CS/vision/grouping/segbench/
count = 0;
ps = 8; % patch size 8 x 8 
fprintf  ('\n%d',length(files));
patches = zeros(ps,ps,length(files)*2400); % 3D array to store patches - size of array is 8 x 8 x N

for i=1:length(files)
    fprintf (' %d',i);
    fname = ['E:\Berkeley\images\train\' files(i).name];
    im = double(rgb2gray(imread(fname)));
    [H,W] = size(im);
    
    for j=1:ps:H-ps+1
        for k=1:ps:W-ps+1
            v = im(j:j+ps-1,k:k+ps-1);
            count = count+1;
            patches(:,:,count) = v;
        end
    end
end

CR = zeros(ps); CC = CR;
for i=1:count
    v = squeeze(patches(:,:,i));
    CC = CC + v*v'; % column-column correlation matrix
    CR = CR + v'*v; % row-row correlation matrix
end
CC = CC/(count-1);
CR = CR/(count-1);

[VR,DR] = eig(CR); % eigenvectors of CR and CC
[VC,DC] = eig(CC);
VR = VR(:,ps:-1:1); % the  eig command of MATLAB sorts the eigenvectors in increasing order of eigenvalues, we want the reverse order
VC = VC(:,ps:-1:1); % the  eig command of MATLAB sorts the eigenvectors in increasing order of eigenvalues, we want the reverse order

VDCT = dctmtx(ps)'; % transpose of the output of dctmtx gives the DCT matrix
fprintf ('\nDCT basis matrix: '); display(VDCT);
fprintf ('\nEigenvectors of row-row correlation matrix: ');  VR
fprintf ('\nEigenvectors of column-column correlation matrix: '); VC

dot_prods = zeros(ps,ps);
for i=1:ps
    for j=1:ps
        dps(i,j) =abs(dot(VDCT(:,i),VR(:,j)));        
    end
end
fprintf ('\nDot products between columns of DCT basis matrix and eigenvectors of row-row correlation matrix: '); display(dps); % show the dot product values
dot_prods = zeros(ps,ps);
for i=1:ps
    for j=1:ps
        dps(i,j) =abs(dot(VDCT(:,i),VC(:,j)));        
    end
end
fprintf ('\nDot products between columns of DCT basis matrix and eigenvectors of column-column correlation matrix: '); display(dps); % show the dot product values

fprintf ('\nRow-row correlation matrix: '); display (CR/CR(1,1)); % see the similarity between this matrix and the one given in the slides under DCT and PCA
fprintf ('\nColumn-column correlation matrix: '); display (CC/CC(1,1)); % see the similarity between this matrix and the one given in the slides under DCT and PCA

figure (1); title ('Basis images: kronecker product of columns of VR and VC');
V = kron(VR,VC);
count  = 0;
for i=1:ps*ps
    count = count+1;
    v = V(:,i); v = (v - min(v(:)))/(max(v(:)) - min(v(:)));    
    subplot(ps,ps,count); imshow(reshape(v,ps,ps));
end

figure (2); title ('DCT Basis images');
V = kron(VDCT,VDCT);
count  = 0;
for i=1:ps*ps
    count = count+1;
    v = V(:,i); if i ~= 1, v = (v - min(v(:)))/(max(v(:)) - min(v(:))); else v = v/max(v(:)); end;
    figure(2); subplot(ps,ps,count); imshow(reshape(v,ps,ps));
end
