Implemented Fisherfaces

.gitignore

@@ -1,4 +1,6 @@
 *.o
 eigenfaces/
-matches/
+fisherfaces/
+matches_eigenfaces/
+matches_fisherfaces/
 main

Eigenfaces.cpp

@@ -24,19 +24,24 @@
 using namespace std;
 using namespace Eigen;
 
-void Eigenfaces::compute(MatrixXf &images) {
-    PCA::compute(images);
+void Eigenfaces::train(const MatrixXf &images) {
+    this->n_pixels = images.rows();
+
+    // Copy values from PCA
+    this->numComponents = PCA::compute(images);
+    this->V = PCA::U; // This contains all Eigenfaces
+
 #ifndef NDEBUG
     cout << "Calculate weights for all images" << endl;
 #endif
-    W_all = project(images);
+    this->W_all = project(images); // Calculate weights
 #ifndef NDEBUG
     cout << "W_all: " << W_all.rows() << " x " << W_all.cols() << endl;
 #endif
-    face_all = reconstructFace(W_all);
+    this->face_all = reconstructFace(W_all);
     //cout << "face_all: " << face_all.rows() << " x " << face_all.cols() << endl;
 }
 
-VectorXf Eigenfaces::euclideanDist(const VectorXf &W) {
-    return ((W_all.colwise() - W)/n_pixels).colwise().norm()/sqrt(K); // Measure euclidean distance between weights
+MatrixXf Eigenfaces::reconstructFace(const MatrixXf &W) {
+    return (V*W).colwise() + PCA::mu; // Reconstruct face
 }

Eigenfaces.h

@@ -15,38 +15,25 @@
  e-mail   :  lauszus@gmail.com
 */
 
-#ifndef _eigenfaces_h_
-#define _eigenfaces_h_
+#ifndef __eigenfaces_h__
+#define __eigenfaces_h__
 
 #include <Eigen/Dense> // http://eigen.tuxfamily.org
 
+#include "Facebase.h"
 #include "PCA.h"
 
 using namespace Eigen;
 
-class Eigenfaces : public PCA {
+class Eigenfaces : public Facebase, public PCA {
 public:
-    Eigenfaces() : PCA() {};
+    void train(const MatrixXf &images);
 
-    void compute(MatrixXf &images);
-    VectorXf euclideanDist(const VectorXf &W);
-
-    template<typename VectorType>
-    vector<size_t> sortIndexes(const VectorType &v) {
-        // Based on: http://stackoverflow.com/a/12399290/2175837
-        // Initialize original index locations
-        vector<size_t> idx(v.size());
-        iota(idx.begin(), idx.end(), 0);
-
-        // Sort indexes based on comparing values in v
-        sort(idx.begin(), idx.end(),
-            [&v](size_t i1, size_t i2) { return v[i1] < v[i2]; });
-
-        return idx;
+    // Facebase implementations
+    MatrixXf project(const MatrixXf &X) {
+        return PCA::project(X); // Simply call PCA implementation
     };
-
-private:
-    MatrixXf W_all, face_all; // Total weights and faces
+    MatrixXf reconstructFace(const MatrixXf &W);
 };
 
 #endif

Facebase.h

@@ -0,0 +1,44 @@
+/* Copyright (C) 2016 Kristian Sloth Lauszus. All rights reserved.
+
+ This software may be distributed and modified under the terms of the GNU
+ General Public License version 2 (GPL2) as published by the Free Software
+ Foundation and appearing in the file GPL2.TXT included in the packaging of
+ this file. Please note that GPL2 Section 2[b] requires that all works based
+ on this software must also be made publicly available under the terms of
+ the GPL2 ("Copyleft").
+
+ Contact information
+ -------------------
+
+ Kristian Sloth Lauszus
+ Web      :  http://www.lauszus.com
+ e-mail   :  lauszus@gmail.com
+*/
+
+#ifndef __facebase_h__
+#define __facebase_h__
+
+#include <Eigen/Dense> // http://eigen.tuxfamily.org
+
+#include "PCA.h"
+
+using namespace Eigen;
+
+class Facebase {
+public:
+    virtual MatrixXf project(const MatrixXf &X) = 0;
+    virtual MatrixXf reconstructFace(const MatrixXf &W) = 0;
+
+    VectorXf euclideanDist(const VectorXf &W) {
+        return ((W_all.colwise() - W)/n_pixels).colwise().norm()/sqrt(numComponents); // Measure euclidean distance between weights
+    };
+
+    MatrixXf V; // Eigenvector
+    int32_t numComponents; // Number of components
+
+protected:
+    MatrixXf W_all, face_all; // Total weights and faces
+    size_t n_pixels;
+};
+
+#endif

Fisherfaces.cpp

@@ -0,0 +1,54 @@
+/* Copyright (C) 2016 Kristian Sloth Lauszus. All rights reserved.
+
+ This software may be distributed and modified under the terms of the GNU
+ General Public License version 2 (GPL2) as published by the Free Software
+ Foundation and appearing in the file GPL2.TXT included in the packaging of
+ this file. Please note that GPL2 Section 2[b] requires that all works based
+ on this software must also be made publicly available under the terms of
+ the GPL2 ("Copyleft").
+
+ Contact information
+ -------------------
+
+ Kristian Sloth Lauszus
+ Web      :  http://www.lauszus.com
+ e-mail   :  lauszus@gmail.com
+*/
+
+#include <iostream>
+
+#include <Eigen/Dense> // http://eigen.tuxfamily.org
+
+#include "Fisherfaces.h"
+
+using namespace std;
+using namespace Eigen;
+
+void Fisherfaces::train(const MatrixXf &images, const VectorXi &classes) {
+    this->n_pixels = images.rows();
+    size_t n_images = classes.size(); // Get number of images
+    int c = classes.maxCoeff(); // Calculate the number of classes, assuming that labels start at 1 and are incremented by 1
+
+    PCA::compute(images, n_images - c);
+    MatrixXf W_pca = PCA::project(images); // Project images onto subspace
+    this->numComponents = LDA::compute(W_pca, classes, c - 1); // Copy number of components from LDA
+
+#ifndef NDEBUG
+    cout << "Calculate weights for all images" << endl;
+#endif
+    this->V = PCA::U*LDA::U; // Calculate Fisherfaces
+    this->W_all = project(images); // Calculate weights
+#ifndef NDEBUG
+    cout << "W_all: " << W_all.rows() << " x " << W_all.cols() << endl;
+#endif
+    this->face_all = reconstructFace(W_all);
+    //cout << "face_all: " << face_all.rows() << " x " << face_all.cols() << endl;
+}
+
+MatrixXf Fisherfaces::project(const MatrixXf &X) {
+    return V.transpose()*X; // Project X onto Fisherface subspace
+}
+
+MatrixXf Fisherfaces::reconstructFace(const MatrixXf &W) {
+    return V*W; // Reconstruct face
+}

Fisherfaces.h

@@ -0,0 +1,38 @@
+/* Copyright (C) 2016 Kristian Sloth Lauszus. All rights reserved.
+
+ This software may be distributed and modified under the terms of the GNU
+ General Public License version 2 (GPL2) as published by the Free Software
+ Foundation and appearing in the file GPL2.TXT included in the packaging of
+ this file. Please note that GPL2 Section 2[b] requires that all works based
+ on this software must also be made publicly available under the terms of
+ the GPL2 ("Copyleft").
+
+ Contact information
+ -------------------
+
+ Kristian Sloth Lauszus
+ Web      :  http://www.lauszus.com
+ e-mail   :  lauszus@gmail.com
+*/
+
+#ifndef __fisherfaces_h__
+#define __fisherfaces_h__
+
+#include <Eigen/Dense> // http://eigen.tuxfamily.org
+
+#include "Facebase.h"
+#include "PCA.h"
+#include "LDA.h"
+
+using namespace Eigen;
+
+class Fisherfaces : public Facebase, public PCA, public LDA {
+public:
+    void train(const MatrixXf &images, const VectorXi &classes);
+
+    // Facebase implementations
+    MatrixXf project(const MatrixXf &X);
+    MatrixXf reconstructFace(const MatrixXf &W);
+};
+
+#endif

LDA.cpp

@@ -0,0 +1,95 @@
+/* Copyright (C) 2016 Kristian Sloth Lauszus. All rights reserved.
+
+ This software may be distributed and modified under the terms of the GNU
+ General Public License version 2 (GPL2) as published by the Free Software
+ Foundation and appearing in the file GPL2.TXT included in the packaging of
+ this file. Please note that GPL2 Section 2[b] requires that all works based
+ on this software must also be made publicly available under the terms of
+ the GPL2 ("Copyleft").
+
+ Contact information
+ -------------------
+
+ Kristian Sloth Lauszus
+ Web      :  http://www.lauszus.com
+ e-mail   :  lauszus@gmail.com
+*/
+
+// Inspired by: https://github.com/opencv/opencv/blob/master/modules/core/src/lda.cpp,
+// https://github.com/bytefish/facerec/blob/master/py/facerec/feature.py,
+// and https://github.com/bytefish/facerec/blob/master/m/models/lda.m
+
+#include <iostream>
+#include <algorithm>
+
+#include <Eigen/Dense> // http://eigen.tuxfamily.org
+#include <Eigen/Eigenvalues>
+#include <RedSVD/RedSVD-h> // https://github.com/ntessore/redsvd-h
+
+#include "LDA.h"
+#include "Tools.h"
+
+using namespace std;
+using namespace Eigen;
+
+// See: http://eigen.tuxfamily.org/dox/structEigen_1_1IOFormat.html
+static IOFormat OctaveFmt(StreamPrecision, 0, ", ", ";\n", "", "", "[", "]");
+
+int32_t LDA::compute(const MatrixXf &X, const VectorXi &classes, int32_t numComponents) {
+#ifndef NDEBUG
+    cout << "Computing LDA" << endl;
+#endif // NDEBUG
+    const size_t N = X.cols();
+    assert((size_t)classes.size() == N); // Make sure that there is a class for each column
+    const size_t dim = X.rows();
+    const int c = classes.maxCoeff(); // Calculate the number of classes, assuming that labels start at 1 and are incremented by 1
+
+#ifndef NDEBUG
+    cout << "N: " << N << endl;
+    cout << "dim: " << dim << endl;
+    cout << "c: " << c << endl;
+#endif // NDEBUG
+
+    numComponents = min(c - 1, numComponents); // Make sure number of components is valid
+#ifndef NDEBUG
+    cout << "numComponents: " << numComponents << endl;
+#endif // NDEBUG
+
+    VectorXf mu = X.rowwise().mean(); // Calculate the mean along each row
+
+    vector<size_t> idx = sortIndexes(classes); // Get indices corresponding to the different classes sorted by the class number
+    //for (auto i : idx) cout << "classes(" << i << "): " << classes(i) << endl;
+
+    MatrixXf Sw = MatrixXf::Zero(dim, dim);
+    MatrixXf Sb = MatrixXf::Zero(dim, dim);
+    // cout << X.format(OctaveFmt) << endl;
+    for (int i = 0; i < c; i++) {
+        int class_count = std::count(classes.data(), classes.data() + classes.size(), i + 1); // Get the number of time that specific class occurs
+        // cout << "class_count: " << class_count << endl;
+
+        MatrixXf X_class(dim, class_count);
+        size_t index = 0;
+        for (auto j : idx) { // TODO: Optimize this, as the indices are already sorted
+            //cout << "classes(" << j << "): " << classes(j) << endl;
+            if (classes(j) == i + 1)
+                X_class.block(0, index++, dim, 1) = X.block(0, j, dim, 1); // Copy data from original image, each pixel at a time based on the indices
+        }
+        // cout << X_class.format(OctaveFmt) << endl;
+
+        VectorXf mu_class = X_class.rowwise().mean(); // Calculate the mean along each row
+        X_class = X_class.colwise() - mu_class; // Subtract means from all columns, thus centering the data
+        VectorXf mu_vector = mu_class - mu;
+
+        Sw += X_class*X_class.transpose(); // Calculate within-class scatter
+        Sb += N*mu_vector*mu_vector.transpose(); // Calculate between-class scatter
+    }
+
+    EigenSolver<MatrixXf> es(Sw.inverse()*Sb); // Solves eigenvalues and eigenvectors of the matrix
+    // cout << es.eigenvalues().head(numComponents).transpose().format(OctaveFmt) << endl;
+    // for (size_t i = 0; i < numComponents; i++)
+    //     cout << es.eigenvectors().col(i).format(OctaveFmt) << endl;
+    U = es.eigenvectors().real().block(0, 0, dim, numComponents); // Extract eigenvectors
+    // cout << U.block(0, 0, dim, 10).format(OctaveFmt) << endl;
+
+    return numComponents;
+}

LDA.h

@@ -0,0 +1,40 @@
+/* Copyright (C) 2016 Kristian Sloth Lauszus. All rights reserved.
+
+ This software may be distributed and modified under the terms of the GNU
+ General Public License version 2 (GPL2) as published by the Free Software
+ Foundation and appearing in the file GPL2.TXT included in the packaging of
+ this file. Please note that GPL2 Section 2[b] requires that all works based
+ on this software must also be made publicly available under the terms of
+ the GPL2 ("Copyleft").
+
+ Contact information
+ -------------------
+
+ Kristian Sloth Lauszus
+ Web      :  http://www.lauszus.com
+ e-mail   :  lauszus@gmail.com
+*/
+
+#ifndef __lda_h__
+#define __lda_h__
+
+#include <Eigen/Dense> // http://eigen.tuxfamily.org
+
+using namespace Eigen;
+
+class LDA {
+public:
+    /**
+     * Computes the Eigenvectors of X using LDA.
+     * @param  X             Data matrix.
+     * @param  classes       Class vector. Must start at 1 and increment by one.
+     * @param  numComponents Number of components used for the analysis.
+     * @return               Returns the number of components used.
+     */
+    int32_t compute(const MatrixXf &X, const VectorXi &classes, int32_t numComponents);
+
+protected:
+    MatrixXf U; // Eigenvectors
+};
+
+#endif

Makefile

@@ -1,8 +1,8 @@
 CC = gcc
 CXX = g++
 
-CFLAGS = -Wall -O3 -std=gnu11 -Wfatal-errors -Wunused-parameter -Wextra
-CXXFLAGS = -Wall -O3 -std=gnu++11 -Wfatal-errors -Wunused-parameter -Wextra
+CFLAGS = -Wall -O3 -std=gnu11 -Wfatal-errors -Wunused-parameter -Wextra -openmp
+CXXFLAGS = -Wall -O3 -std=gnu++11 -Wfatal-errors -Wunused-parameter -Wextra -openmp
 
 # Include Eigen library
 CXXFLAGS += `pkg-config --cflags eigen3`