Download the code from Matlab

Contents

Image compression by Reduce Basis Decomposition

Author: Yanlai Chen

Version: 1.0

Date: 03/20/2015

thispic = imread('Figures/lena.png');
compratio = 10;
figshrink = 1;

m = size(thispic,1);
n = size(thispic,2);
num = floor(min(m,n)/compratio);
fprintf('Shrinking from (%d, %d) to (%d, %d) and (%d, %d)\n',m,n,m,num,num,n);
datad3 = size(thispic,3);

Shrinking from (512, 512) to (512, 51) and (51, 512)

Displaying the original picture as is

figure('Position',[1,1,1+n/figshrink,1+m/figshrink])
if(datad3 == 1)
    imagesc(thispic);colormap(gray);drawnow;
else
    image(thispic);
    title('Original image');drawnow;
end
axis off

Handling it by RBD on each component

bases = zeros(m,num,datad3);
TransMat = zeros(num,n,datad3);
tic
for i=1:datad3
    data = double(thispic(:,:,i));
    [bases(:,:,i), TransMat(:,:,i)] = RBD(data, 1e-20,num);
end
t_rb = toc;

Recovering the compressed image from its RBD bases and displaying it

B10 = zeros(m,n,datad3);
for j=1:datad3
    B10(:,:,j) = bases(:,:,j)*TransMat(:,:,j);
end
ScaleAndSave(B10,m,n,datad3,'By RB component-wise compression',figshrink);

Handling it by SVDS

U = zeros(m,num,datad3);
S = zeros(num,num,datad3);
V = zeros(n,num,datad3);

U = zeros(m,num,datad3);
S = zeros(num,num,datad3);
V = zeros(n,num,datad3);

tic
for i=1:datad3
    data = double(thispic(:,:,i));
    [U(:,:,i), S(:,:,i), V(:,:,i)] = svds(data,num);
end
t_svds = toc;

Handling it by SVD

U = zeros(m,m,datad3);
S = zeros(m,n,datad3);
V = zeros(n,n,datad3);
tic
for i=1:datad3
    data = double(thispic(:,:,i));
    [U(:,:,i), S(:,:,i), V(:,:,i)] = svd(data);
end
t_svd = toc;

Recovering the compressed image from its SVDD bases and displaying it

A10 = zeros(m,n,datad3);
for j=1:datad3
    for i=1:num
        A10(:,:,j) = A10(:,:,j) + S(i,i,j)*U(:,i,j)*V(:,i,j)';
    end
end
ScaleAndSave(A10,m,n,datad3,'By SVD',figshrink);

RBD image compression is faster than SVD

fprintf('Time spend on (svd, svds, rbd) relative to svd is (%e, %e, %e)\n',t_svd/t_svd, t_svds/t_svd, t_rb/t_svd);

Time spend on (svd, svds, rbd) relative to svd is (1.000000e+00, 3.907092e+00, 2.636303e-01)

Face recognition by Reduced Basis Decomposition.

We test our algorithm on the UMIST data set available at Roweis' website http://www.cs.nyu.edu/ roweis/data.html

load('Figures/umist_cropped.mat');
totalp = 20;

It contains 575 cropped facial images (20 people under different poses) Here is a snapshot of the database

umist_snap = imread('Figures/umist_snapshot.png');
figure('Position',[1,1,1+size(umist_snap,2),1+size(umist_snap,1)])
imagesc(umist_snap);colormap(gray);drawnow
axis off

Preliminary setup

trainnum: number of views randomly take per subject for training.

totalrun: number of rounds of formations of training and test sets

basnum: number of total eigen faces taken from the whole training set

avgrate: the average success rate when totalrun rounds are done

trainnum = 10;
totalrun = 2;
avgrate = [];
for basnum = 5:10:100
    thisrate = 0;
    for runnum = 1: totalrun

        trainpic2d = [];
        tp_ind = zeros(totalp,trainnum);
        testpic2d = [];

Formation of training set and test set

We first lineup all the pictures and then randomly form the training set and test set.

        for i=1:totalp
            pic_pi = facedat{i};
            num_pi = size(pic_pi,3);
            [nr,nc] = size(pic_pi(:,:,1));

            pic2d = zeros(nr*nc,num_pi);
            for j=1:num_pi
                pic2d(:,j) = reshape(pic_pi(:,:,j),nr*nc,1);
            end
            pic2ds{i} = pic2d;

            tmp = 0;
            while (tmp < trainnum)
                tn = randi(num_pi);
                if(size(find(tp_ind(i,1:tmp) == tn),2) == 0)
                    tp_ind(i,tmp+1) = tn;
                    tmp = tmp+1;
                    trainpic2d = [trainpic2d double(reshape(pic_pi(:,:,tn),nr*nc,1))];
                end
            end
        end

Method number 1: A simple PCA approach to find the eigen faces

        tic
        PCA_FR(trainpic2d);
        t_pca = toc;

Method number 2: RBD to find the (basnum) eigen faces

        tic
        totalmean = mean(trainpic2d,2);
        trainpic2d = trainpic2d - totalmean*ones(1,size(trainpic2d,2));
        [bases, TransMat] = RBD(trainpic2d, 1e-20,basnum,1,zeros(nr*nc,1));
        t_rbd = toc;

The comparison in efficiency

        fprintf('PCA time over RBD time when %d eigen faces are used is %f\n',basnum,t_pca/t_rbd);

PCA time over RBD time when 5 eigen faces are used is 2.608177

PCA time over RBD time when 15 eigen faces are used is 1.893146

PCA time over RBD time when 25 eigen faces are used is 1.483418

PCA time over RBD time when 35 eigen faces are used is 1.071684

PCA time over RBD time when 45 eigen faces are used is 0.845395

PCA time over RBD time when 55 eigen faces are used is 0.702518

PCA time over RBD time when 65 eigen faces are used is 0.559466

PCA time over RBD time when 75 eigen faces are used is 0.450852

PCA time over RBD time when 85 eigen faces are used is 0.372011

PCA time over RBD time when 95 eigen faces are used is 0.321043

Testing the eigen faces and calculate the success rate

totaltest: the total number of tests you want to perform

        matched = 0;
        totaltest = 100;
        for i=1:totaltest
            pnum = randi(totalp);
            %    picnum = tp_ind(pnum,randi(trainnum)); % test those in the
            %    training set, success rate should be 100%

            % test those that are not in the training set
            intraining = 1;
            while(intraining == 1)
                picnum = randi(size(pic2ds{pnum},2));
                if(size(find(tp_ind(pnum,:) == picnum),2) == 0)
                    intraining = 0;
                end
            end

            thispic = pic2ds{pnum}(:,picnum) - totalmean;

            all_err = 1e14;
            bestall = 0;
            coe = bases'*thispic;
            for j=1:size(TransMat,2)
                jth_err = norm(coe - TransMat(:,j));
                if(jth_err <all_err)
                    all_err = jth_err;
                    bestall = j;
                end
            end
            bestall = ceil(bestall/(trainnum+1e-10));

            if(bestall == pnum)
                matched = matched+1;
            end
        end
       thisrate = thisrate + matched/totaltest;

    end

    thisrate = thisrate/totalrun;
    fprintf('The average success rate using %d eigen faces for %d runs of recognizing %d images is %f\n',basnum,totalrun,totaltest,thisrate)
    avgrate = [avgrate thisrate];
end

The average success rate using 5 eigen faces for 2 runs of recognizing 100 images is 0.795000
The average success rate using 15 eigen faces for 2 runs of recognizing 100 images is 0.930000
The average success rate using 25 eigen faces for 2 runs of recognizing 100 images is 0.940000
The average success rate using 35 eigen faces for 2 runs of recognizing 100 images is 0.960000
The average success rate using 45 eigen faces for 2 runs of recognizing 100 images is 0.910000
The average success rate using 55 eigen faces for 2 runs of recognizing 100 images is 0.950000
The average success rate using 65 eigen faces for 2 runs of recognizing 100 images is 0.915000
The average success rate using 75 eigen faces for 2 runs of recognizing 100 images is 0.905000
The average success rate using 85 eigen faces for 2 runs of recognizing 100 images is 0.950000
The average success rate using 95 eigen faces for 2 runs of recognizing 100 images is 0.980000

Plotting the error rates

figure()
plot(5:10:100,100 - avgrate*100,'-r','LineWidth',2,'MarkerEdgeColor','k','MarkerFaceColor','g','MarkerSize',10);
xlabel('Number of eigen faces','fontsize',15);
ylabel('Classification error rate (%)','fontsize',15);

Published with MATLAB® R2014a

Contents

function []= PCA_FR(AllImgs)

The training part of a typical face recognition algorithm using PCA

Author: Yanlai Chen

Date: 03/22/2015

% Preprocessing of the images
m = mean(AllImgs,2);
imgcount = size(AllImgs,2);

SubM = [];
for i=1 : imgcount
    temp = double(AllImgs(:,i)) - m;
    SubM = [SubM temp];
end

The Eigen solve

that is usually needed by all PCA-based algorithms V : eigenvector matrix D : eigenvalue matrix

L = SubM' * SubM;
[V,D]=eig(L);

% Kaiser's rule on finding how many Principal Components to take
KEV = [];
for i = 1 : size(V,2)
    if( D(i,i) > 1 )
        KEV = [KEV V(:,i)];
    end
end

% building the EigenFaces and preparing for the test phase
eigenfaces = SubM * KEV;

ProjImg = [ ];
for i = 1 : size(eigenfaces,2)
    temp = eigenfaces' * SubM(:,i);
    ProjImg = [ProjImg temp];
end

end

Published with MATLAB® R2014a

Contents

function [bases, TransMat] = RBD(data, tol, col, NormMat, StartingVec)

Reduced basis decomposition

This function approximate a matrix (data) by a product of two matrices (bases and TransMat) with the column space of the first matrix approximating that of the input matrix.

There are multiple ways to call it as you can specify from 1 to 5 inputs

RBD(data), RBD(data,tol), RBD(data,tol,col) etc

% Inputs:
%
% Provide only the first three unless you understand what you are doing
%
% * data => the data matrix you want to decompose
% * tol  => the accuracy you desire of your decomposition
% * col  => the number of columns you can afford to have in the compressed
%           matrix
% * NormMat => A SPD matrix if you wish to use a different norm to measure
%           the error
% * StartingVec => The vector you wish to use to start the greedy algorithm
%
% Outputs:
%
% bases*TransMat $\approx$ data
%
% * bases    => the "compressed" matrix (to be more accurate, the group of vectors spanning
% approximately the column space of data)
% * TransMat => the "transformation" matrix
%
% Author: Yanlai Chen
%
% Version: 1.0
%
% Date: 03/20/2015

Error Checking

We provide values to those inputs that are not supplied

if(nargin<1)
    error('RBD needs at least one input, a matrix')
end

if (nargin == 1)
    tol = 1e-6;
    col = min(10,size(data,2));
    TrueE_NoR = 1;
    StartingVec = 0;
elseif(nargin == 2)
    col = size(data,2);
    TrueE_NoR = 1;
    StartingVec = 0;
elseif(nargin == 3)
    TrueE_NoR = 1;
    StartingVec = 0;
elseif(nargin == 4)
    TrueE_NoR = 0;
    StartingVec = 0;
else
    TrueE_NoR = 0;
end

Praparation of the algorithm

nr = size(data,1);
nc = size(data,2);

bases = zeros(nr,col);

TransMat = zeros(col,nc);
xiFlag=zeros(1,col);
xiFlag(1) = randi(nc);
i=1;
CurErr = tol + 1;

% Preparation for efficient error evaluation
ftf = sum(data'.*data',2);
AtAAtXi = zeros(nc,col);

if(TrueE_NoR == 0)
    AtAAt = data'*NormMat;
    XitAAtXi = zeros(col,col);
    tM = data'*NormMat;
% a very fast way to evaluate ftf(j) = data(:,j)'*AAt*data(:,j);
    ftf = sum(tM.*data',2);
end

The RBD greedy algorithm

while (i <= col) && (CurErr > tol)

    if((i==1) && norm(StartingVec) > 1e-6)
        biCand = StartingVec;
    else
        biCand = data(:,xiFlag(i));
    end

Inside: Gram-Schmidt orthonormalization of the current candidate with all

previsouly chosen basis vectors

    for j=1:i-1
        biCand = biCand - (biCand'*bases(:,j))*bases(:,j);
    end
    normi = sqrt(biCand'*biCand);
    if(normi < 1e-7)
        fprintf('Reduced system getting singular - to stop with %d basis functions\n',i-1);
        bases = bases(:,1:i-1);
        TransMat = TransMat(1:i-1,:);
        break
    else
        bases(:,i) = biCand/normi;
    end
    TransMat(i,:) = bases(:,i)'*data;

Inside: With one more basis added, we need to update what allows for the

efficient error evaluation.

    if(TrueE_NoR == 0)
        AtAAtXi(:,i) = AtAAt*bases(:,i);
        XitAAtXi(i,1:i) = bases(:,i)'*NormMat*bases(:,1:i);
        XitAAtXi(1:i,i) = XitAAtXi(i,1:i)';
    else
        AtAAtXi(:,i) = data'*bases(:,i);
    end

Inside: Efficiently go through all the columns to identify where the error would

be the largest if we were to use the current space for compression.

    TMM = TransMat(1:i,:);
    te1 = sum(AtAAtXi(:,1:i).*TMM',2);
    if(TrueE_NoR == 0)
        tM = TMM'*XitAAtXi(1:i,1:i);
    else
        tM = TMM';
    end
    te2 = sum(tM.*TMM',2);
    errord = ftf - 2*te1 + te2;
    [CurErr,TempPos] = max(errord);

% Mark this location for the next round
    if(i < col)
        xiFlag(i+1) = TempPos;
    end
% If the largest error is small enough, we announce and stop.
    if(CurErr <= tol)
        fprintf('Reduced system getting accurate enough - to stop with %d basis functions\n',i);
        bases = bases(:,1:i);
        TransMat = TransMat(1:i,:);
    else
        i = i+1;
    end

end

end

Published with MATLAB® R2014a

Contents

    function [] = ScaleAndSave(pictosave,m,n,datad3,picname,figshrink)

A simple code for scaling and saving a picture

Author: Yanlai Chen

Date: 03/22/2015

    for sj=1:datad3
            if(datad3 == 3)
                jmin = min(min(pictosave(:,:,sj)));
                jmax = max(max(pictosave(:,:,sj)));
                pictosave(:,:,sj) = (pictosave(:,:,sj) - jmin)/(jmax - jmin);
            end
        end

        figure('Position',[1,1,1+n/figshrink,1+m/figshrink])
        if(datad3 == 1)
            imagesc(pictosave);colormap(gray);drawnow;
        else
            image(pictosave);
            title(picname);
            drawnow;
        end
        axis off

    end

Published with MATLAB® R2014a