Principal Component Analysis matlab code

The Principal Component Analysis (PCA) is one of the most successful techniques that have been used in image recognition and compression. PCA is a statistical method under the broad title of factor analysis. The purpose of PCA is to reduce the large dimensionality of the data space (observed variables) to the smaller intrinsic dimensionality of feature space (independent variables), which are needed to describe the data economically. This is the case when there is a strong correlation between observed variables.

The jobs which PCA can do are prediction, redundancy removal, feature extraction, data compression, etc. Because PCA is a classical technique which can do something in the linear domain, applications having linear models are suitable, such as signal processing, image processing, system and control theory, communications, etc



WHEN PUBLISHING A PAPER AS A RESULT OF RESEARCH CONDUCTED BY USING THIS CODE
OR ANY PART OF IT, MAKE A REFERENCE TO THE FOLLOWING PAPER:
Delac K., Grgic M., Grgic S., Independent Comparative Study of PCA, ICA, and LDA
on the FERET Data Set, International Journal of Imaging Systems and Technology,
Vol. 15, Issue 5, 2006, pp. 252-260

=====================================================================================

pca.m

-------------------------------------------------------------------------------------

function pca (path, trainList, subDim)

if nargin < 3
subDim = dim - 1;
end;

disp(' ')

load listAll;

% Constants
numIm = 3816;


% Memory allocation for DATA matrix
fprintf('Creating DATA matrix\n')
tmp = imread ( [path char(listAll(1)) '.pgm'] );
[m, n] = size (tmp); % image size - used later also!!!
DATA = uint8 (zeros(m*n, numIm)); % Memory allocated
clear str tmp;

% Creating DATA matrix
for i = 1 : numIm
im = imread ( [path char(listAll(i)) '.pgm'] );
DATA(:, i) = reshape (im, m*n, 1);
end;
save DATA DATA;
clear im;

% Creating training images space
fprintf('Creating training images space\n')
dim = length (trainList);
imSpace = zeros (m*n, dim);
for i = 1 : dim
index = strmatch (trainList(i), listAll);
imSpace(:, i) = DATA(:, index);
end;
save imSpace imSpace;
clear DATA;

% Calculating mean face from training images
fprintf('Zero mean\n')
psi = mean(double(imSpace'))';
save psi psi;

% Zero mean
zeroMeanSpace = zeros(size(imSpace));
for i = 1 : dim
zeroMeanSpace(:, i) = double(imSpace(:, i)) - psi;
end;
save zeroMeanSpace zeroMeanSpace;
clear imSpace;

% PCA
fprintf('PCA\n')
L = zeroMeanSpace' * zeroMeanSpace; % Turk-Pentland trick (part 1)
[eigVecs, eigVals] = eig(L);

diagonal = diag(eigVals);
[diagonal, index] = sort(diagonal);
index = flipud(index);

pcaEigVals = zeros(size(eigVals));
for i = 1 : size(eigVals, 1)
pcaEigVals(i, i) = eigVals(index(i), index(i));
pcaEigVecs(:, i) = eigVecs(:, index(i));
end;

pcaEigVals = diag(pcaEigVals);
pcaEigVals = pcaEigVals / (dim-1);
pcaEigVals = pcaEigVals(1 : subDim); % Retaining only the largest subDim ones

pcaEigVecs = zeroMeanSpace * pcaEigVecs; % Turk-Pentland trick (part 2)

save pcaEigVals pcaEigVals;

% Normalisation to unit length
fprintf('Normalising\n')
for i = 1 : dim
pcaEigVecs(:, i) = pcaEigVecs(:, i) / norm(pcaEigVecs(:, i));
end;

% Dimensionality reduction.
fprintf('Creating lower dimensional subspace\n')
w = pcaEigVecs(:, 1:subDim);
save w w;
clear w;

% Subtract mean face from all images
load DATA;
load psi;
zeroMeanDATA = zeros(size(DATA));
for i = 1 : size(DATA, 2)
zeroMeanDATA(:, i) = double(DATA(:, i)) - psi;
end;
clear psi;
clear DATA;

% Project all images onto a new lower dimensional subspace (w)
fprintf('Projecting all images onto a new lower dimensional subspace\n')
load w;
pcaProj = w' * zeroMeanDATA;
clear w;
clear zeroMeanDATA;
save pcaProj pcaProj;

=====================================================================================

createDistMat.m

-------------------------------------------------------------------------------------

function distMat = createDistMat (proj, metric)

% Memory allocation
distMat = zeros(max(size(proj)));

switch (metric)

case 'L1'
distMat = pdist(proj', 'cityblock');
case 'L2'
distMat = pdist(proj', 'euclidean');
case 'COS'
distMat = pdist(proj', 'cosine');

otherwise
error('%s metric not supported.', metric);

end; % switch (metric) ends here

distMat = squareform(distMat);

=====================================================================================

feret.m

-------------------------------------------------------------------------------------

function [FERET] = feret(distMat, rank)

% Load feretGallery and probe lists
load listAll;
load feretGallery;
load fb;
load fc;
load dup1;
load dup2;

% Constants (number of images in feretGallery and probes)
numGalleryImgs = size (feretGallery, 1);
numFbImgs = size (fb, 1);
numFcImgs = size (fc, 1);
numDup1Imgs = size (dup1, 1);
numDup2Imgs = size (dup2, 1);

% Get the list of positions where feretGallery images are located in the list
% of all images and store it in variable index
feretGalleryIndex = getIndex (feretGallery, listAll);

% Get the list of all the probe images
fbIndex = getIndex (fb, listAll);
fcIndex = getIndex (fc, listAll);
dup1Index = getIndex (dup1, listAll);
dup2Index = getIndex (dup2, listAll);

% Calculate ranks for the CMS curve
% The results are stores in a structure
FERET.fb = getResults (distMat, feretGallery, feretGalleryIndex, fb, fbIndex, numFbImgs, rank);
FERET.fc = getResults (distMat, feretGallery, feretGalleryIndex, fc, fcIndex, numFcImgs, rank);
FERET.dup1 = getResults (distMat, feretGallery, feretGalleryIndex, dup1, dup1Index, numDup1Imgs, rank);
FERET.dup2 = getResults (distMat, feretGallery, feretGalleryIndex, dup2, dup2Index, numDup2Imgs, rank);

%**************************************************************************
% Statistics (CMS curve, total number of probe images)
%**************************************************************************

function [RESULTS] = getResults (distMat, feretGallery, feretGalleryIndex, probe, probeIndex, numProbeImgs, rank)

for j = 1 : rank
for i = 1 : numProbeImgs

position = probeIndex(i);
currentRow = distMat(position,:);
reduced = currentRow(1, feretGalleryIndex);

[Y, I] = sort(reduced);
inx = I(j);

% Determine if G=H based on the first 5 characters of the filename
G = char(feretGallery(inx));
H = char(probe(i));
correct(i) = strncmp(G, H, 5);
% if G=H correct(i)=1, else 0

end;

% Rank 1 result (percentage)
if j == 1
RESULTS.rank1 = (sum(correct)/numProbeImgs)*100;
end;

% Percentage of correctly recognized images for a given probe set
RESULTS.perc(j) = (sum(correct)/numProbeImgs)*100;

% CMS curve
RESULTS.cms(j) = sum(RESULTS.perc(1:j));

end;

% Total number of probe images
RESULTS.numProbeImgs = numProbeImgs;
%**************************************************************************

%**************************************************************************
% Find positions of probe or feretGallery images in the list of all images
%**************************************************************************

function [index] = getIndex (sub, all)

num = size (sub, 1);
for i = 1 : num
index(i) = strmatch(sub(i), all);
end;
%**************************************************************************

Viola-Jones object detection framework

The Viola-Jones object detection framework is the first object detection framework to provide competitive object detection rates in real-time proposed in 2001 by Paul Viola and Michael Jones.

Introduction
Object detection is detecting a specified object class such as cars, faces, plates ext. in a given image or a video sequence. Object detection has many applications in computer based vision such as object tracking, object recognition, and scene surveillance.

The technique relies on the use of simple Haar-like features that are evaluated quickly through the use of a new image representation. Based on the concept of an “Integral Image” it generates a large set of features and uses the boosting algorithm AdaBoost to reduce the over-complete set and the introduction of a degenerative tree of the boosted classifiers provides for robust and fast interferences. The detector is applied in a scanning fashion and used on gray-scale images, the scanned window that is applied can also be scaled, as well as the features evaluated.

In the technique only simple rectangular (Haar-like) features are used, reminiscent to Haar basis functions. These features are equivalent to intensity difference readings and are quite easy to compute. There are three feature types used with varying numbers of sub-rectangles, two, two rectangles, one three and one four rectangle feature types. Using rectangular features instead of the pixels in an image provides a number of benefits, namely a sort of a ad-hoc domain knowledge is implied as well as a speed increase over pixel based systems. The calculation of the features is facilitated with the use of an “integral image”. With the introduction of a integral image Viola and Jones are able to calculate in one pass of the sample image, and is one of the keys to the speed of the system. An integral image is similar to a “summed are table”, used in computer graphics but its use is applied in pixel area evaluation.

It was outlined that the implementation of a system that used such features would provide a feature set that was far too large, hence the feature set must be only restricted to a small number of critical features. This is done with the use of boosting algorithm, AdaBoost. Interference is enhanced with the use of AdaBoost where a small set of features is selected from a large set, and in doing so a strong hypothesis is formed, in this case resulting in a strong classifier. Simply having a reduced set of features was not enough to reduce the vast amounts of computation in a detector task, since it is naturally a probabilistic one, hence Viola and Jones proposed the use of degenerative tree of classifiers.

Described by Viola and Jones as a degenerative tree, and sometimes referred to as a decision stump, its use also speeds the detection process. A degenerative tree is the daisy chaining of general to specific classifiers, whereby the first few classifiers are general enough to discount an image sub window and so on the time of further observations by the more specific classifiers down the chain, this can save a large degree of computation.

Integral Image
In order to be successful a face detection algorithm must possess two key features, accuracy and speed. There is generally a trade-off between the two. Through the use of a new image representation, termed integral images, Viola and Jones describe a means for fast feature evaluation, and this proves to be an effective means to speed up the classification task of the system.

Integral images are easy to understand, they are constructed by simply taking the sum of the luminance values above and to the left of a pixel in an image. Viola and Jones make note of the fact that the integral image is effectively the double integral of the sample image, first along the rows then along the columns. Integral images are equivalent to summed-area tables, yet their use is not texture mapping, being so, their implementation us quite well documented.

1 1 1
1 1 1
1 1 1


1 2 3
2 4 6
3 6 9
The brilliance in using an integral image to speed up a feature extraction lies in the fact that any rectangle in an image can be calculated from that images integral image, in only four indexes to the integral image. This makes the otherwise exhaustive process of summing luminances quite rapid. In fact the calculation of an images integral image can be calculated in only one pass of the image, and Matlab experiments have shown that a large set of images (12000) can be calculated within less than 2 seconds.

Integral application
Given a rectangle specified as four coordinates A(x1,y1) upper left and D(x4,y4) lower right, evaluating the area of the rectangle is done in four image references to the integral image, this represents a huge performance increase in terms of feature extraction.

Sum of grey rectangle = D - (B + C) + A

Since both rectangle B and C include rectangle A the sum of A has to be added to the calculation.


Source:
http://www.codeproject.com/Articles/85113/Efficient-Face-Detection-Algorithm-using-Viola-Jon.aspx

Robust Face Detection in C/C++ (Haar-like features)

Best solution might be Haar-like features.

Viola and Jones adapted the idea of using Haar wavelets and developed the so called Haar-like features. A Haar-like feature considers adjacent rectangular regions at a specific location in a detection window, sums up the pixel intensities in these regions and calculates the difference between them. This difference is then used to categorize subsections of an image. For example, let us say we have an image database with human faces. It is a common observation that among all faces the region of the eyes is darker than the region of the cheeks. Therefore a common haar feature for face detection is a set of two adjacent rectangles that lie above the eye and the cheek region. The position of these rectangles is defined relative to a detection window that acts like a bounding box to the target object (the face in this case).

source: http://en.wikipedia.org/wiki/Haar-like_features

Different Algorithims of Face Recognition

PCA
Derived from Karhunen-Loeve's transformation. Given an s-dimensional vector representation of each face in a training set of images, Principal Component Analysis (PCA) tends to find a t-dimensional subspace whose basis vectors correspond to the maximum variance direction in the original image space. This new subspace is normally lower dimensional .If the image elements are considered as random variables, the PCA basis vectors are defined as eigenvectors of the scatter matrix.

ICA
Independent Component Analysis (ICA) minimizes both second-order and higher-order dependencies in the input data and attempts to find the basis along which the data (when projected onto them) are - statistically independent . Bartlett et al. provided two architectures of ICA for face recognition task: Architecture I - statistically independent basis images, and Architecture II - factorial code representation.

LDA
Linear Discriminant Analysis (LDA) finds the vectors in the underlying space that best discriminate among classes. For all samples of all classes the between-class scatter matrix SB and the within-class scatter matrix SW are defined. The goal is to maximize SB while minimizing SW, in other words, maximize the ratio detSB/detSW . This ratio is maximized when the column vectors of the projection matrix are the eigenvectors of (SW^-1 × SB).

EP
Aa eigenspace-based adaptive approach that searches for the best set of projection axes in order to maximize a fitness function, measuring at the same time the classification accuracy and generalization ability of the system. Because the dimension of the solution space of this problem is too big, it is solved using a specific kind of genetic algorithm called Evolutionary Pursuit (EP).

EBGM
Elastic Bunch Graph Matching (EBGM). All human faces share a similar topological structure. Faces are represented as graphs, with nodes positioned at fiducial points. (exes, nose...) and edges labeled with 2-D distance vectors. Each node contains a set of 40 complex Gabor wavelet coefficients at different scales and orientations (phase, amplitude). They are called "jets". Recognition is based on labeled graphs. A labeled graph is a set of nodes connected by edges, nodes are labeled with jets, edges are labeled with distances.

Kernel Methods
The face manifold in subspace need not be linear. Kernel methods are a generalization of linear methods. Direct non-linear manifold schemes are explored to learn this non-linear manifold.

Trace Transform
The Trace transform, a generalization of the Radon transform, is a new tool for image processing which can be used for recognizing objects under transformations, e.g. rotation, translation and scaling. To produce the Trace transform one computes a functional along tracing lines of an image. Different Trace transforms can be produced from an image using different trace functionals.

AAM
An Active Appearance Model (AAM) is an integrated statistical model which combines a model of shape variation with a model of the appearance variations in a shape-normalized frame. An AAM contains a statistical model if the shape and gray-level appearance of the object of interest which can generalize to almost any valid example. Matching to an image involves finding model parameters which minimize the difference between the image and a synthesized model example projected into the image.

3-D Morphable Model
Human face is a surface lying in the 3-D space intrinsically. Therefore the 3-D model should be better for representing faces, especially to handle facial variations, such as pose, illumination etc. Blantz et al. proposed a method based on a 3-D morphable face model that encodes shape and texture in terms of model parameters, and algorithm that recovers these parameters from a single image of a face.

3-D Face Recognition
The main novelty of this approach is the ability to compare surfaces independent of natural deformations resulting from facial expressions. First, the range image and the texture of the face are acquired. Next, the range image is preprocessed by removing certain parts such as hair, which can complicate the recognition process. Finally, a canonical form of the facial surface is computed. Such a representation is insensitive to head orientations and facial expressions, thus significantly simplifying the recognition procedure. The recognition itself is performed on the canonical surfaces.

Digital Camera Face Recognition: How It Works?

Face detection technology, available from manufacturers such as Canon, Pentax and Fuji­Film, uses special algorithms to parse the scene while you aim the camera. What's it looking for? The shape of a human face, of course. When it finds one, the camera automatically adjusts both the focus and the exposure to provide the best portrait possible. In the case of FujiFilm's "Image Intelligence" system, for example, a chip inside the camera constantly scans the image in its viewfinder for two eyes, a nose, ears and a chin, making out up to 10 faces at a time before you've hit the shutter. Usually, a scanned face has to cover at least 10 percent of the height of the viewfinder LCD — this requirement prevents the camera from trying to lock on to faces far in the background. The face identification process takes about four one-hundredths of a second. You can turn off the face recognition feature for times when your subject isn't a human. However, according to FujiFilm project manager David Troy, research shows that "70 percent of images taken have a person as the subject." Face recognition can be used for more than just autofocus. Some cameras capture the location of the identified faces within each picture snapped. This lets you zoom in on the faces automatically on the camera's LCD after you take the shot — useful when checking if grandma's eyes were closed or open. Sometime soon, face detection may even give way to facial identification, discerning one subject from an-other. For instance, the camera could retain an image tagged "Mom" in its memory. Later, the camera would automatically recognize each subsequent picture of your mother and add the "Mom" tag to it. As for the question about animal faces, you're not alone in wondering how the new cameras handle our friends from other species. "You'd be surprised how often we get that question!" says FujiFilm spokeswoman Katherine Keane. "Unfortunately, FujiFilm's face detection technology can only detect human faces." So you'll just have to focus on Fido's friendly mug yourself.

Laplace edge detection with C++ code

The 5x5 Laplacian used is a convoluted mask to approximate the second derivative, unlike the Sobel method which approximates the gradient. And instead of 2 3x3 Sobel masks, one for the x and y direction, Laplace uses 1 5x5 mask for the 2nd derivative in both the x and y directions. However, because these masks are approximating a second derivative measurement on the image, they are very sensitive to noise, as can be seen by comparing edgeSob.bmp to edgeLap.bmp. The Laplace mask and code are shown below:

for(Y=0; Y<=(originalImage.rows-1); Y++) {
for(X=0; X<=(originalImage.cols-1); X++) {
SUM = 0;

/* image boundaries */
if(Y==0 Y==1 Y==originalImage.rows-2 Y==originalImage.rows-1)
SUM = 0;
else if(X==0 X==1 X==originalImage.cols-2 X==originalImage.cols-1)
SUM = 0;

/* Convolution starts here */
else {
for(I=-2; I<=2; I++) {
for(J=-2; J<=2; J++) {
SUM = SUM + (int)( (*(originalImage.data + X + I +
(Y + J)*originalImage.cols)) * MASK[I+2][J+2]);

}
}
}
if(SUM>255) SUM=255;
if(SUM<0) SUM=0;

*(edgeImage.data + X + Y*originalImage.cols) = 255 - (unsigned char)(SUM);
fwrite((edgeImage.data + X + Y*originalImage.cols),sizeof(char),1,bmpOutput);
}
}