diff --git a/app/br/br.cpp b/app/br/br.cpp
index 9e79f37..2e2293a 100644
--- a/app/br/br.cpp
+++ b/app/br/br.cpp
@@ -163,8 +163,8 @@ public:
                 check((parc >= 2) && (parc <= 6), "Incorrect parameter count for 'evalDetection'.");
                 br_eval_detection(parv[0], parv[1], parc >= 3 ? parv[2] : "", parc >= 4 ? atoi(parv[3]) : 0, parc >= 5 ? atoi(parv[4]) : 0, parc == 6 ? atoi(parv[5]) : 0);
             } else if (!strcmp(fun, "evalLandmarking")) {
-                check((parc >= 2) && (parc <= 5), "Incorrect parameter count for 'evalLandmarking'.");
-                br_eval_landmarking(parv[0], parv[1], parc >= 3 ? parv[2] : "", parc >= 4 ? atoi(parv[3]) : 0, parc >= 5 ? atoi(parv[4]) : 1);
+                check((parc >= 2) && (parc <= 7), "Incorrect parameter count for 'evalLandmarking'.");
+                br_eval_landmarking(parv[0], parv[1], parc >= 3 ? parv[2] : "", parc >= 4 ? atoi(parv[3]) : 0, parc >= 5 ? atoi(parv[4]) : 1,  parc >= 6 ? atoi(parv[5]) : 0, parc >= 7 ? atoi(parv[6]) : 5);
             } else if (!strcmp(fun, "evalRegression")) {
                 check(parc >= 2 && parc <= 4, "Incorrect parameter count for 'evalRegression'.");
                 br_eval_regression(parv[0], parv[1], parc >= 3 ? parv[2] : "", parc >= 4 ? parv[3] : "");
@@ -264,7 +264,7 @@ private:
                "-evalClassification <predicted_gallery> <truth_gallery> <predicted property name> <ground truth proprty name>\n"
                "-evalClustering <clusters> <gallery>\n"
                "-evalDetection <predicted_gallery> <truth_gallery> [{csv}] [{normalize}] [{minSize}]\n"
-               "-evalLandmarking <predicted_gallery> <truth_gallery> [{csv} [<normalization_index_a> <normalization_index_b>]]\n"
+               "-evalLandmarking <predicted_gallery> <truth_gallery> [{csv} [<normalization_index_a> <normalization_index_b>] [sample_index] [total_examples]]\n"
                "-evalRegression <predicted_gallery> <truth_gallery> <predicted property name> <ground truth property name>\n"
                "-assertEval <simmat> <mask> <accuracy>\n"
                "-plotDetection <file> ... <file> {destination}\n"
diff --git a/openbr/core/common.cpp b/openbr/core/common.cpp
index d6c445a..c3802e4 100644
--- a/openbr/core/common.cpp
+++ b/openbr/core/common.cpp
@@ -16,21 +16,32 @@
 
 #include "common.h"
 #include <QMutex>
+#include <RandomLib/Random.hpp>
 
 using namespace std;
 
+static RandomLib::Random g_rand;
+static QMutex rngLock;
+
 /**** GLOBAL ****/
 void Common::seedRNG() {
-    static QMutex seedControl;
-    QMutexLocker lock(&seedControl);
+    QMutexLocker lock(&rngLock);
 
     static bool seeded = false;
     if (!seeded) {
         srand(0); // We seed with 0 instead of time(NULL) to have reproducible randomness
         seeded = true;
+        g_rand.Reseed(0);
     }
 }
 
+double Common::randN()
+{
+    QMutexLocker lock(&rngLock);
+
+    return g_rand.FloatN();
+}
+
 QList<int> Common::RandSample(int n, int max, int min, bool unique)
 {
     QList<int> samples; samples.reserve(n);
diff --git a/openbr/core/common.h b/openbr/core/common.h
index 1a9be82..843dc96 100644
--- a/openbr/core/common.h
+++ b/openbr/core/common.h
@@ -220,6 +220,9 @@ double KernelDensityEstimation(const V<T> &vals, double x, double h)
     return y / (vals.size() * h);
 }
 
+// Return a random number, uniformly distributed over 0,1
+double randN();
+
 /*!
  * \brief Returns a vector of n integers sampled in the range <min, max].
  *
@@ -236,19 +239,14 @@ QList<int> RandSample(int n, const QSet<int> &values, bool unique = false);
 template <typename T>
 QList<int> RandSample(int n, const QList<T> &weights, bool unique = false)
 {
-    static bool seeded = false;
-    if (!seeded) {
-        srand(time(NULL));
-        seeded = true;
-    }
-
     QList<T> cdf = CumSum(weights);
     for (int i=0; i<cdf.size(); i++) // Normalize cdf
         cdf[i] = cdf[i] / cdf.last();
 
     QList<int> samples; samples.reserve(n);
     while (samples.size() < n) {
-        T r = (T)rand() / (T)RAND_MAX;
+        T r = randN();
+
         for (int j=0; j<weights.size(); j++) {
             if ((r >= cdf[j]) && (r <= cdf[j+1])) {
                 if (!unique || !samples.contains(j))
diff --git a/openbr/core/core.cpp b/openbr/core/core.cpp
index 4d4bb01..a6f68bf 100644
--- a/openbr/core/core.cpp
+++ b/openbr/core/core.cpp
@@ -110,9 +110,11 @@ struct AlgorithmCore
 
     void store(const QString &model) const
     {
-        // Create stream
-        QByteArray data;
-        QDataStream out(&data, QFile::WriteOnly);
+        QtUtils::BlockCompression compressedWrite;
+        QFile outFile(model);
+        compressedWrite.setBasis(&outFile);
+        QDataStream out(&compressedWrite);
+        compressedWrite.open(QFile::WriteOnly);
 
         // Serialize algorithm to stream
         transform->serialize(out);
@@ -131,18 +133,16 @@ struct AlgorithmCore
         if (mode == TransformCompare)
             comparison->serialize(out);
 
-        // Compress and save to file
-        QtUtils::writeFile(model, data, -1);
+        compressedWrite.close();
     }
 
     void load(const QString &model)
     {
-        // Load from file and decompress
-        QByteArray data;
-        QtUtils::readFile(model, data, true);
-
-        // Create stream
-        QDataStream in(&data, QFile::ReadOnly);
+        QtUtils::BlockCompression compressedRead;
+        QFile inFile(model);
+        compressedRead.setBasis(&inFile);
+        QDataStream in(&compressedRead);
+        compressedRead.open(QFile::ReadOnly);
 
         // Load algorithm
         transform = QSharedPointer<Transform>(Transform::deserialize(in));
diff --git a/openbr/core/eigenutils.cpp b/openbr/core/eigenutils.cpp
index c701d44..b1ff605 100644
--- a/openbr/core/eigenutils.cpp
+++ b/openbr/core/eigenutils.cpp
@@ -4,80 +4,15 @@
 using namespace Eigen;
 using namespace cv;
 
-//Helper function to quickly write eigen matrix to disk. Not efficient.
-void writeEigen(MatrixXf X, QString filename) {
-    Mat m(X.rows(),X.cols(),CV_32FC1);
-    for (int i = 0; i < X.rows(); i++) {
-        for (int j = 0; j < X.cols(); j++) {
-            m.at<float>(i,j) = X(i,j);
-        }
-    }
-    QScopedPointer<br::Format> format(br::Factory<br::Format>::make(filename));
-    format->write(br::Template(m));
-}
-
-void writeEigen(MatrixXd X, QString filename) {
-    Mat m(X.rows(),X.cols(),CV_32FC1);
-    for (int i = 0; i < X.rows(); i++) {
-        for (int j = 0; j < X.cols(); j++) {
-            m.at<float>(i,j) = (float)X(i,j);
-        }
-    }
-    QScopedPointer<br::Format> format(br::Factory<br::Format>::make(filename));
-    format->write(br::Template(m));
-}
-
-void writeEigen(VectorXd X, QString filename) {
-    Mat m(X.size(),1,CV_32FC1);
-    for (int i = 0; i < X.rows(); i++) {
-        m.at<float>(i,0) = (float)X(i);
-    }
-    QScopedPointer<br::Format> format(br::Factory<br::Format>::make(filename));
-    format->write(br::Template(m));
-}
-
-void writeEigen(VectorXf X, QString filename) {
-    Mat m(X.size(),1,CV_32FC1);
-    for (int i = 0; i < X.rows(); i++) {
-        m.at<float>(i,0) = X(i);
-    }
-    QScopedPointer<br::Format> format(br::Factory<br::Format>::make(filename));
-    format->write(br::Template(m));
-}
-
-void printEigen(Eigen::MatrixXd X) {
-    for (int i = 0; i < X.rows(); i++) {
-        QString str;
-        for (int j = 0; j < X.cols(); j++) {
-            str.append(QString::number(X(i,j)) + " ");
-        }
-        qDebug() << str;
-    }
-}
-void printEigen(Eigen::MatrixXf X) {
-    for (int i = 0; i < X.rows(); i++) {
-        QString str;
-        for (int j = 0; j < X.cols(); j++) {
-            str.append(QString::number(X(i,j)) + " ");
-        }
-        qDebug() << str;
-    }
-}
-
-void printSize(Eigen::MatrixXf X) {
+void EigenUtils::printSize(Eigen::MatrixXf X) {
     qDebug() << "Rows=" << X.rows() << "\tCols=" << X.cols();
 }
 
-float eigMean(const Eigen::MatrixXf& x) {
-    return x.array().sum() / (x.rows() * x.cols());
+float EigenUtils::stddev(const Eigen::MatrixXf& x) {
+    return sqrt((x.array() - x.mean()).pow(2).sum() / (x.cols() * x.rows()));
 }
 
-float eigStd(const Eigen::MatrixXf& x) {
-    float mean = eigMean(x);
-    return sqrt((x.array() - mean).pow(2).sum() / (x.cols() * x.rows()));
-}
-
-MatrixXf removeRowCol(const MatrixXf X, int row, int col) {
+MatrixXf EigenUtils::removeRowCol(const MatrixXf X, int row, int col) {
     MatrixXf Y(X.rows() - 1,X.cols() - 1);
 
     for (int i1 = 0, i2 = 0; i1 < X.rows(); i1++) {
@@ -96,7 +31,7 @@ MatrixXf removeRowCol(const MatrixXf X, int row, int col) {
     return Y;
 }
 
-MatrixXf pointsToMatrix(const QList<QPointF> points, bool isAffine) {
+MatrixXf EigenUtils::pointsToMatrix(const QList<QPointF> points, bool isAffine) {
     MatrixXf P(points.size(), isAffine ? 3 : 2);
     for (int i = 0; i < points.size(); i++) {
         P(i, 0) = points[i].x();
@@ -107,7 +42,7 @@ MatrixXf pointsToMatrix(const QList<QPointF> points, bool isAffine) {
     return P;
 }
 
-QList<QPointF> matrixToPoints(const Eigen::MatrixXf P) {
+QList<QPointF> EigenUtils::matrixToPoints(const Eigen::MatrixXf P) {
     QList<QPointF> points;
     for (int i = 0; i < P.rows(); i++)
         points.append(QPointF(P(i, 0), P(i, 1)));
@@ -115,7 +50,7 @@ QList<QPointF> matrixToPoints(const Eigen::MatrixXf P) {
 }
 
 //Converts x y points in a single vector to two column matrix
-Eigen::MatrixXf vectorToMatrix(const Eigen::MatrixXf vector) {
+Eigen::MatrixXf EigenUtils::vectorToMatrix(const Eigen::MatrixXf vector) {
     int n = vector.rows();
     Eigen::MatrixXf matrix(n / 2, 2);
     for (int i = 0; i < n / 2; i++) {
@@ -126,7 +61,7 @@ Eigen::MatrixXf vectorToMatrix(const Eigen::MatrixXf vector) {
     return matrix;
 }
 
-Eigen::MatrixXf matrixToVector(const Eigen::MatrixXf matrix) {
+Eigen::MatrixXf EigenUtils::matrixToVector(const Eigen::MatrixXf matrix) {
     int n2 = matrix.rows();
     Eigen::MatrixXf vector(n2 * 2, 1);
     for (int i = 0; i < n2; i++) {
@@ -136,12 +71,3 @@ Eigen::MatrixXf matrixToVector(const Eigen::MatrixXf matrix) {
     }
     return vector;
 }
-
-Eigen::MatrixXf toEigen(const Mat m) {
-    if (m.type() != CV_32F)
-        qFatal("Mat to Eigen Converstation only supports CV_32F");
-
-    Eigen::MatrixXf data(m.rows, m.cols);
-    return Eigen::Map<const Eigen::MatrixXf>(m.ptr<float>(), m.rows, m.cols);
-}
-
diff --git a/openbr/core/eigenutils.h b/openbr/core/eigenutils.h
index 73a72e0..fae0511 100644
--- a/openbr/core/eigenutils.h
+++ b/openbr/core/eigenutils.h
@@ -18,32 +18,74 @@
 #define EIGENUTILS_H
 
 #include <QDataStream>
+#include <QDebug>
 #include <Eigen/Core>
 #include <assert.h>
 
+#include <opencv2/core/eigen.hpp>
 #include <opencv2/core/core.hpp>
 
-void writeEigen(Eigen::MatrixXf X, QString filename);
-void writeEigen(Eigen::MatrixXd X, QString filename);
-void writeEigen(Eigen::VectorXd X, QString filename);
-void writeEigen(Eigen::VectorXf X, QString filename);
-void printEigen(Eigen::MatrixXd X);
-void printEigen(Eigen::MatrixXf X);
-void printSize(Eigen::MatrixXf X);
+#include "openbr/core/qtutils.h"
 
-//Converts x y points in a single vector to two column matrix
-Eigen::MatrixXf vectorToMatrix(const Eigen::MatrixXf vector);
-Eigen::MatrixXf matrixToVector(const Eigen::MatrixXf matrix);
+namespace EigenUtils
+{
+    template<typename _Scalar, int _Rows, int _Cols, int _Options, int _MaxRows, int _MaxCols>
+    QString matrixToString(const Eigen::Matrix< _Scalar, _Rows, _Cols, _Options, _MaxRows, _MaxCols > &mat)
+    {
+        QString result;
+        if (mat.rows() > 1) result += "{ ";
+        for (int r=0; r<mat.rows(); r++) {
+            if ((mat.rows() > 1) && (r > 0)) result += "  ";
+            if (mat.cols() > 1) result += "[";
+            for (int c=0; c<mat.cols(); c++) {
+                result += QString::number(mat(r, c));
+                if (c < mat.cols() - 1) result += ", ";
+            }
+            if (mat.cols() > 1) result += "]";
+            if (r < mat.rows()-1) result += "\n";
+        }
+        if (mat.rows() > 1) result += " }";
+        return result;
+    }
+
+    template<typename _Scalar, int _Rows, int _Cols, int _Options, int _MaxRows, int _MaxCols>
+    void writeMatrix(const Eigen::Matrix< _Scalar, _Rows, _Cols, _Options, _MaxRows, _MaxCols > &mat, const QString &filename)
+    {
+        int r = mat.rows();
+        int c = mat.cols();
+
+        _Scalar *data = new _Scalar[r*c];
+        for (int i=0; i<r; i++)
+            for (int j=0; j<c; j++)
+                data[i*c+j] = mat(i, j);
+        int bytes = r*c*sizeof(_Scalar);
+        QByteArray byteArray((const char*)data,bytes);
+        QtUtils::writeFile(filename,byteArray);
+    }
+
+    void printSize(Eigen::MatrixXf X);
 
-//Remove row and column from the matrix:
-Eigen::MatrixXf removeRowCol(const Eigen::MatrixXf X, int row, int col);
+    // Converts x y points in a single vector to two column matrix
+    Eigen::MatrixXf vectorToMatrix(const Eigen::MatrixXf vector);
+    Eigen::MatrixXf matrixToVector(const Eigen::MatrixXf matrix);
 
-//Convert a point list into a matrix:
-Eigen::MatrixXf pointsToMatrix(const QList<QPointF> points, bool isAffine=false);
-QList<QPointF> matrixToPoints(const Eigen::MatrixXf P);
+    // Remove row and column from the matrix:
+    Eigen::MatrixXf removeRowCol(const Eigen::MatrixXf X, int row, int col);
 
-//Convert cv::Mat to Eigen
-Eigen::MatrixXf toEigen(const cv::Mat m);
+    // Convert a point list into a matrix:
+    Eigen::MatrixXf pointsToMatrix(const QList<QPointF> points, bool isAffine=false);
+    QList<QPointF> matrixToPoints(const Eigen::MatrixXf P);
+
+    // Compute the element-wise standard deviation
+    float stddev(const Eigen::MatrixXf& x);
+}
+
+template<typename _Scalar, int _Rows, int _Cols, int _Options, int _MaxRows, int _MaxCols>
+inline QDebug operator<<(QDebug dbg, const Eigen::Matrix< _Scalar, _Rows, _Cols, _Options, _MaxRows, _MaxCols > &mat)
+{
+    dbg.nospace() << EigenUtils::matrixToString(mat);
+    return dbg.space();
+}
 
 template<typename _Scalar, int _Rows, int _Cols, int _Options, int _MaxRows, int _MaxCols>
 inline QDataStream &operator<<(QDataStream &stream, const Eigen::Matrix< _Scalar, _Rows, _Cols, _Options, _MaxRows, _MaxCols > &mat)
@@ -83,56 +125,4 @@ inline QDataStream &operator>>(QDataStream &stream, Eigen::Matrix< _Scalar, _Row
     return stream;
 }
 
-/*Compute the mean of the each column (dim == 1) or row (dim == 2)
-  of the matrix*/
-template<typename T>
-Eigen::MatrixBase<T> eigMean(const Eigen::MatrixBase<T>& x,int dim)
-{
-    if (dim == 1) {
-        Eigen::MatrixBase<T> y(1,x.cols());
-        for (int i = 0; i < x.cols(); i++)
-            y(i) = x.col(i).sum() / x.rows();
-        return y;
-    } else if (dim == 2) {
-        Eigen::MatrixBase<T> y(x.rows(),1);
-        for (int i = 0; i < x.rows(); i++)
-            y(i) = x.row(i).sum() / x.cols();
-        return y;
-    }
-    qFatal("A matrix can only have two dimensions");
-}
-
-/*Compute the element-wise mean*/
-float eigMean(const Eigen::MatrixXf& x);
-/*Compute the element-wise mean*/
-float eigStd(const Eigen::MatrixXf& x);
-
-/*Compute the std dev of the each column (dim == 1) or row (dim == 2)
-  of the matrix*/
-template<typename T>
-Eigen::MatrixBase<T> eigStd(const Eigen::MatrixBase<T>& x,int dim)
-{
-    Eigen::MatrixBase<T> mean = eigMean(x, dim);
-    if (dim == 1) {
-        Eigen::MatrixBase<T> y(1,x.cols());
-        for (int i = 0; i < x.cols(); i++) {
-            T value = 0;
-            for (int j = 0; j < x.rows(); j++)
-                value += pow(y(j, i) - mean(i), 2);
-            y(i) = sqrt(value / (x.rows() - 1));
-        }
-        return y;
-    } else if (dim == 2) {
-        Eigen::MatrixBase<T> y(x.rows(),1);
-        for (int i = 0; i < x.rows(); i++) {
-            T value = 0;
-            for (int j = 0; j < x.cols(); j++)
-                value += pow(y(i, j) - mean(j), 2);
-            y(i) = sqrt(value / (x.cols() - 1));
-        }
-        return y;
-    }
-    qFatal("A matrix can only have two dimensions");
-}
-
 #endif // EIGENUTILS_H
diff --git a/openbr/core/eval.cpp b/openbr/core/eval.cpp
index ef5bc1c..006b322 100755
--- a/openbr/core/eval.cpp
+++ b/openbr/core/eval.cpp
@@ -18,6 +18,7 @@
 #include "eval.h"
 #include "openbr/core/common.h"
 #include "openbr/core/qtutils.h"
+#include "openbr/core/opencvutils.h"
 #include <QMapIterator>
 
 using namespace cv;
@@ -1068,16 +1069,25 @@ float EvalDetection(const QString &predictedGallery, const QString &truthGallery
     return averageOverlap;
 }
 
-float EvalLandmarking(const QString &predictedGallery, const QString &truthGallery, const QString &csv, int normalizationIndexA, int normalizationIndexB)
+static void projectAndWrite(Transform *t, const Template &src, const QString &filePath)
+{
+    Template dst;
+    t->project(src,dst);
+    OpenCVUtils::saveImage(dst.m(),filePath);
+}
+
+float EvalLandmarking(const QString &predictedGallery, const QString &truthGallery, const QString &csv, int normalizationIndexA, int normalizationIndexB, int sampleIndex, int totalExamples)
 {
     qDebug("Evaluating landmarking of %s against %s", qPrintable(predictedGallery), qPrintable(truthGallery));
-    const TemplateList predicted(TemplateList::fromGallery(predictedGallery));
-    const TemplateList truth(TemplateList::fromGallery(truthGallery));
-    const QStringList predictedNames = File::get<QString>(predicted, "name");
-    const QStringList truthNames = File::get<QString>(truth, "name");
+    TemplateList predicted(TemplateList::fromGallery(predictedGallery));
+    TemplateList truth(TemplateList::fromGallery(truthGallery));
+    QStringList predictedNames = File::get<QString>(predicted, "name");
+    QStringList truthNames = File::get<QString>(truth, "name");
 
     int skipped = 0;
     QList< QList<float> > pointErrors;
+    QList<float> imageErrors;
+    QList<float> normalizedLengths;
     for (int i=0; i<predicted.size(); i++) {
         const QString &predictedName = predictedNames[i];
         const int truthIndex = truthNames.indexOf(predictedName);
@@ -1085,39 +1095,94 @@ float EvalLandmarking(const QString &predictedGallery, const QString &truthGalle
         const QList<QPointF> predictedPoints = predicted[i].file.points();
         const QList<QPointF> truthPoints = truth[truthIndex].file.points();
         if (predictedPoints.size() != truthPoints.size() || truthPoints.contains(QPointF(-1,-1))) {
-            skipped++;
+            predicted.removeAt(i);
+            predictedNames.removeAt(i);
+            truth.removeAt(i);
+            truthNames.removeAt(i);
+            i--; skipped++;
             continue;
         }
+
         while (pointErrors.size() < predictedPoints.size())
             pointErrors.append(QList<float>());
+
+        // Want to know error for every image.
+
         if (normalizationIndexA >= truthPoints.size()) qFatal("Normalization index A is out of range.");
         if (normalizationIndexB >= truthPoints.size()) qFatal("Normalization index B is out of range.");
         const float normalizedLength = QtUtils::euclideanLength(truthPoints[normalizationIndexB] - truthPoints[normalizationIndexA]);
-        for (int j=0; j<predictedPoints.size(); j++)
-            pointErrors[j].append(QtUtils::euclideanLength(predictedPoints[j] - truthPoints[j])/normalizedLength);
+        normalizedLengths.append(normalizedLength);
+        float totalError = 0;
+        for (int j=0; j<predictedPoints.size(); j++) {
+            float error = QtUtils::euclideanLength(predictedPoints[j] - truthPoints[j])/normalizedLength;
+            totalError += error;
+            pointErrors[j].append(error);
+        }
+        imageErrors.append(totalError/predictedPoints.size());
     }
-    qDebug() << "Skipped " << skipped << " files due to point size mismatch.";
+
+    qDebug() << "Skipped" << skipped << "files due to point size mismatch.";
 
     QList<float> averagePointErrors; averagePointErrors.reserve(pointErrors.size());
-    for (int i=0; i<pointErrors.size(); i++) {
-        std::sort(pointErrors[i].begin(), pointErrors[i].end());
-        averagePointErrors.append(Common::Mean(pointErrors[i]));
-    }
-    const float averagePointError = Common::Mean(averagePointErrors);
 
     QStringList lines;
     lines.append("Plot,X,Y");
+
+    QtUtils::touchDir(QDir("landmarking_examples_truth"));
+    QtUtils::touchDir(QDir("landmarking_examples_predicted"));
+
+    // Example
+    {
+        QScopedPointer<Transform> t(Transform::make("Open+Draw(verbose,rects=false,location=false)",NULL));
+
+        QString filePath = "landmarking_examples_truth/"+truth[sampleIndex].file.fileName();
+        projectAndWrite(t.data(), truth[sampleIndex],filePath);
+        lines.append("Sample,"+filePath+","+QString::number(truth[sampleIndex].file.points().size()));
+    }
+
+    // Get best and worst performing examples
+    QList< QPair<float,int> > exampleIndices = Common::Sort(imageErrors,true);
+
+    QScopedPointer<Transform> t(Transform::make("Open+Draw(rects=false)",NULL));
+
+    for (int i=0; i<totalExamples; i++) {
+        QString filePath = "landmarking_examples_truth/"+truth[exampleIndices[i].second].file.fileName();
+        projectAndWrite(t.data(), truth[exampleIndices[i].second],filePath);
+        lines.append("EXT,"+filePath+","+QString::number(exampleIndices[i].first));
+
+        filePath = "landmarking_examples_predicted/"+predicted[exampleIndices[i].second].file.fileName();
+        projectAndWrite(t.data(), predicted[exampleIndices[i].second],filePath);
+        lines.append("EXP,"+filePath+","+QString::number(exampleIndices[i].first));
+    }
+
+    for (int i=exampleIndices.size()-1; i>exampleIndices.size()-totalExamples-1; i--) {
+        QString filePath = "landmarking_examples_truth/"+truth[exampleIndices[i].second].file.fileName();
+        projectAndWrite(t.data(), truth[exampleIndices[i].second],filePath);
+        lines.append("EXT,"+filePath+","+QString::number(exampleIndices[i].first));
+
+        filePath = "landmarking_examples_predicted/"+predicted[exampleIndices[i].second].file.fileName();
+        projectAndWrite(t.data(), predicted[exampleIndices[i].second],filePath);
+        lines.append("EXP,"+filePath+","+QString::number(exampleIndices[i].first));
+    }
+
     for (int i=0; i<pointErrors.size(); i++) {
+        std::sort(pointErrors[i].begin(), pointErrors[i].end());
+        averagePointErrors.append(Common::Mean(pointErrors[i]));
         const QList<float> &pointError = pointErrors[i];
         const int keep = qMin(Max_Points, pointError.size());
         for (int j=0; j<keep; j++)
             lines.append(QString("Box,%1,%2").arg(QString::number(i), QString::number(pointError[j*(pointError.size()-1)/(keep-1)])));
     }
 
+    const float averagePointError = Common::Mean(averagePointErrors);
+
     lines.append(QString("AvgError,0,%1").arg(averagePointError));
+    lines.append(QString("NormLength,0,%1").arg(Common::Mean(normalizedLengths)));
 
     QtUtils::writeFile(csv, lines);
-    qDebug("Average Error: %.3f", averagePointError);
+
+    qDebug("Average Error for all Points: %.3f", averagePointError);
+
     return averagePointError;
 }
 
diff --git a/openbr/core/eval.h b/openbr/core/eval.h
index bd83c98..a150e92 100644
--- a/openbr/core/eval.h
+++ b/openbr/core/eval.h
@@ -31,7 +31,7 @@ namespace br
 
     void EvalClassification(const QString &predictedGallery, const QString &truthGallery, QString predictedProperty = "", QString truthProperty = "");
     float EvalDetection(const QString &predictedGallery, const QString &truthGallery, const QString &csv = "", bool normalize = false, int minSize = 0, int maxSize = 0); // Return average overlap
-    float EvalLandmarking(const QString &predictedGallery, const QString &truthGallery, const QString &csv = "", int normalizationIndexA = 0, int normalizationIndexB = 1); // Return average error
+    float EvalLandmarking(const QString &predictedGallery, const QString &truthGallery, const QString &csv = "", int normalizationIndexA = 0, int normalizationIndexB = 1, int sampleIndex = 0, int totalExamples = 5); // Return average error
     void EvalRegression(const QString &predictedGallery, const QString &truthGallery, QString predictedProperty = "", QString truthProperty = "");
 }
 
diff --git a/openbr/core/plot.cpp b/openbr/core/plot.cpp
index 288382d..26a39ee 100644
--- a/openbr/core/plot.cpp
+++ b/openbr/core/plot.cpp
@@ -482,16 +482,115 @@ bool PlotLandmarking(const QStringList &files, const File &destination, bool sho
     qDebug("Plotting %d landmarking file(s) to %s", files.size(), qPrintable(destination));
     RPlot p(files, destination, false);
 
-    p.file.write("# Split data into individual plots\n"
-                 "plot_index = which(names(data)==\"Plot\")\n"
-                 "Box <- data[grep(\"Box\",data$Plot),-c(1)]\n"
-                 "rm(data)\n"
-                 "\n");
-
-    p.file.write(qPrintable(QString("ggplot(Box, aes(Y,%1%2))").arg(p.major.size > 1 ? QString(", colour=%1").arg(p.major.header) : QString(), p.minor.size > 1 ? QString(", linetype=%1").arg(p.minor.header) : QString()) +
+    p.file.write(qPrintable(QString("# Split data into individual plots\n"
+                                    "plot_index = which(names(data)==\"Plot\")\n"
+                                    "Box <- data[grep(\"Box\",data$Plot),-c(1)]\n"
+                                    "Box$X <- factor(Box$X, levels = Box$X, ordered = TRUE)\n"
+                                    "Sample <- data[grep(\"Sample\",data$Plot),-c(1)]\n"
+                                    "Sample$X <- as.character(Sample$X)\n"
+                                    "EXT <- data[grep(\"EXT\",data$Plot),-c(1)]\n"
+                                    "EXT$X <- as.character(EXT$X)\n"
+                                    "EXP <- data[grep(\"EXP\",data$Plot),-c(1)]\n"
+                                    "EXP$X <- as.character(EXP$X)\n"
+                                    "NormLength <- data[grep(\"NormLength\",data$Plot),-c(1)]\n"
+                                    "rm(data)\n"
+                                    "\n")));
+
+    p.file.write(qPrintable(QString("summarySE <- function(data=NULL, measurevar, groupvars=NULL, na.rm=FALSE, conf.interval=.95, .drop=TRUE) {\n\t"
+                                    "require(plyr)\n\n\tlength2 <- function (x, na.rm=FALSE) {\n\t\tif (na.rm) sum(!is.na(x))\n\t\telse       length(x)"
+                                    "\n\t}\n\n\tdatac <- ddply(data, groupvars, .drop=.drop, .fun = function(xx, col) {\n\t\t"
+                                    "c(N=length2(xx[[col]], na.rm=na.rm), mean=mean(xx[[col]], na.rm=na.rm), sd=sd(xx[[col]], na.rm=na.rm))\n\t\t},"
+                                    "\n\t\tmeasurevar\n\t)\n\n\tdatac <- rename(datac, c(\"mean\" = measurevar))\n\tdatac$se <- datac$sd / sqrt(datac$N)"
+                                    "\n\tciMult <- qt(conf.interval/2 + .5, datac$N-1)\n\tdatac$ci <- datac$se * ciMult\n\n\treturn(datac)\n}\n")));
+
+
+    p.file.write(qPrintable(QString("\nreadData <- function(data) {\n\texamples <- list()\n"
+                                    "\tfor (i in 1:nrow(data)) {\n"
+                                    "\t\tpath <- data[i,1]\n"
+                                    "\t\tvalue <- data[i,2]\n"
+                                    "\t\tfile <- unlist(strsplit(path, \"[.]\"))[1]\n"
+                                    "\t\text <- unlist(strsplit(path, \"[.]\"))[2]\n"
+                                    "\t\tif (ext == \"jpg\" || ext == \"JPEG\" || ext == \"jpeg\" || ext == \"JPG\") {\n"
+                                    "\t\t\timg <- readJPEG(path)\n"
+                                    "\t\t} else if (ext == \"PNG\" || ext == \"png\") {\n"
+                                    "\t\t\timg <- readPNG(path)\n"
+                                    "\t\t} else if (ext == \"TIFF\" || ext == \"tiff\" || ext == \"TIF\" || ext == \"tif\") { \n"
+                                    "\t\t\timg <- readTIFF(path)\n"
+                                    "}else {\n"
+                                    "\t\t\tnext\n"
+                                    "\t\t}\n"
+                                    "\t\texample <- list(file = file, value = value, image = img)\n"
+                                    "\t\texamples[[i]] <- example\n"
+                                    "\t}\n"
+                                    "\treturn(examples)\n"
+                                    "}\n")));
+
+    p.file.write(qPrintable(QString("\nlibrary(jpeg)\n"
+                                    "library(png)\n"
+                                    "library(grid)\n"
+                                    "multiplot <- function(..., plotlist=NULL, cols) {\n"
+                                    "\trequire(grid)\n"
+                                    "\t# Make a list from the ... arguments and plotlist\n"
+                                    "\tplots <- c(list(...), plotlist)\n"
+                                    "\tnumPlots = length(plots)\n"
+                                    "\t# Make the panel\n"
+                                    "\tplotCols = cols\n"
+                                    "\tplotRows = ceiling(numPlots/plotCols)\n"
+                                    "\t# Set up the page\n"
+                                    "\tgrid.newpage()\n"
+                                    "\tpushViewport(viewport(layout = grid.layout(plotRows, plotCols)))\n"
+                                    "\tvplayout <- function(x, y)\n"
+                                    "\tviewport(layout.pos.row = x, layout.pos.col = y)\n"
+                                    "\t# Make each plot, in the correct location\n"
+                                    "\tfor (i in 1:numPlots) {\n"
+                                    "\t\tcurRow = ceiling(i/plotCols)\n"
+                                    "\t\tcurCol = (i-1) %% plotCols + 1\n"
+                                    "\t\tprint(plots[[i]], vp = vplayout(curRow, curCol))\n"
+                                    "\t}\n"
+                                    "}\n")));
+
+    p.file.write(qPrintable(QString("\nplotImage <- function(image, title=NULL, label=NULL) { \n"
+                                    "\tp <- qplot(1:10, 1:10, geom=\"blank\") + annotation_custom(rasterGrob(image$image), xmin=-Inf, xmax=Inf, ymin=-Inf, ymax=Inf) + theme(axis.line=element_blank(), axis.title.y=element_blank(), axis.text.x=element_blank(), axis.text.y=element_blank(), line=element_blank(), axis.ticks=element_blank(), panel.background=element_blank()) + labs(title=title) + xlab(label)\n"
+                                    "\treturn(p)"
+                                    "}\n")));
+
+    p.file.write(qPrintable(QString("\nsample <- readData(Sample) \n"
+                                    "rows <- sample[[1]]$value\n"
+                                    "algs <- unique(Box$%1)\n"
+                                    "algs <- algs[!duplicated(algs)]\n"
+                                    "print(plotImage(sample[[1]],\"Sample Landmarks\",sprintf(\"Total Landmarks: %s\",sample[[1]]$value))) \n"
+                                    "if (nrow(EXT) != 0 && nrow(EXP)) {\n"
+                                    "\tfor (j in 1:length(algs)) {\n"
+                                    "\ttruthSample <- readData(EXT[EXT$. == algs[[j]],])\n"
+                                    "\tpredictedSample <- readData(EXP[EXP$. == algs[[j]],])\n"
+                                    "\t\tfor (i in 1:length(predictedSample)) {\n"
+                                    "\t\t\tmultiplot(plotImage(predictedSample[[i]],sprintf(\"%s\\nPredicted Landmarks\",algs[[j]]),sprintf(\"Average Landmark Error: %.3f\",predictedSample[[i]]$value)),plotImage(truthSample[[i]],\"Ground Truth\\nLandmarks\",\"\"),cols=2)\n"
+                                    "\t\t}\n"
+                                    "\t}\n"
+                                    "}\n").arg(p.major.size > 1 ? p.major.header : (p.minor.header.isEmpty() ? p.major.header : p.minor.header))));
+
+    p.file.write(qPrintable(QString("\n"
+                 "# Code to format error table\n"
+                 "StatBox <- summarySE(Box, measurevar=\"Y\", groupvars=c(\"%1\",\"X\"))\n"
+                 "OverallStatBox <- summarySE(Box, measurevar=\"Y\", groupvars=c(\"%1\"))\n"
+                 "mat <- matrix(paste(as.character(round(StatBox$Y, 3)), round(StatBox$ci, 3), sep=\" \\u00b1 \"),nrow=rows,ncol=length(algs),byrow=FALSE)\n"
+                 "mat <- rbind(mat, paste(as.character(round(OverallStatBox$Y, 3)), round(OverallStatBox$ci, 3), sep=\" \\u00b1 \"))\n"
+                 "mat <- rbind(mat, as.character(round(NormLength$Y, 3)))\n"
+                 "colnames(mat) <- algs\n"
+                 "rownames(mat) <- c(seq(0,rows-1),\"Aggregate\",\"Average IPD\")\n"
+                 "ETable <- as.table(mat)\n").arg(p.major.size > 1 ? p.major.header : (p.minor.header.isEmpty() ? p.major.header : p.minor.header))));
+
+    p.file.write(qPrintable(QString("\n"
+                       "print(textplot(ETable))\n"
+                       "print(title(\"Landmarking Error Rates\"))\n")));
+
+    p.file.write(qPrintable(QString("ggplot(Box, aes(Y,%1%2))").arg(p.major.size > 1 ? QString(", colour=%1").arg(p.major.header) : QString(),
+                                                                    p.minor.size > 1 ? QString(", linetype=%1").arg(p.minor.header) : QString()) +
                             QString(" + annotation_logticks(sides=\"b\") + stat_ecdf() + scale_x_log10(\"Normalized Error\", breaks=c(0.001,0.01,0.1,1,10)) + scale_y_continuous(\"Cumulative Density\", label=percent) + theme_minimal()\n\n")));
+
     p.file.write(qPrintable(QString("ggplot(Box, aes(factor(X), Y%1%2))").arg(p.major.size > 1 ? QString(", colour=%1").arg(p.major.header) : QString(), p.minor.size > 1 ? QString(", linetype=%1").arg(p.minor.header) : QString()) +
-                            QString("+ annotation_logticks(sides=\"l\") + geom_boxplot(alpha=0.5) + geom_jitter(size=1, alpha=0.5) + scale_x_discrete(\"Landmark\") + scale_y_log10(\"Normalized Error\", breaks=c(0.01,0.1,1,10)) + theme_minimal()\n\n")));
+                            QString("+ annotation_logticks(sides=\"l\") + geom_boxplot(alpha=0.5) + geom_jitter(size=1, alpha=0.5) + scale_x_discrete(\"Landmark\") + scale_y_log10(\"Normalized Error\", breaks=c(0.001,0.01,0.1,1,10)) + theme_minimal()\n\n")));
+
     p.file.write(qPrintable(QString("ggplot(Box, aes(factor(X), Y%1%2))").arg(p.major.size > 1 ? QString(", colour=%1").arg(p.major.header) : QString(), p.minor.size > 1 ? QString(", linetype=%1").arg(p.minor.header) : QString()) +
                             QString("+ annotation_logticks(sides=\"l\") + geom_violin(alpha=0.5) + scale_x_discrete(\"Landmark\") + scale_y_log10(\"Normalized Error\", breaks=c(0.001,0.01,0.1,1,10))\n\n")));
 
diff --git a/openbr/core/qtutils.cpp b/openbr/core/qtutils.cpp
index 0f189b0..fdfd004 100644
--- a/openbr/core/qtutils.cpp
+++ b/openbr/core/qtutils.cpp
@@ -500,6 +500,131 @@ QString getAbsolutePath(const QString &filename)
     return QFileInfo(filename).absoluteFilePath();
 }
 
+BlockCompression::BlockCompression(QIODevice *_basis)
+{
+    blockSize = 100000000;
+    setBasis(_basis);
+}
+
+BlockCompression::BlockCompression() { blockSize = 100000000; };
+
+
+bool BlockCompression::open(QIODevice::OpenMode mode)
+{
+    this->setOpenMode(mode);
+    bool res = basis->open(mode);
+
+    if (!res)
+        return false;
+
+    blockReader.setDevice(basis);
+    blockWriter.setDevice(basis);
+
+    if (mode & QIODevice::WriteOnly) {
+        precompressedBlockWriter = new QBuffer;
+        precompressedBlockWriter->open(QIODevice::ReadWrite);
+    }
+    else if (mode & QIODevice::ReadOnly) {
+        QByteArray compressedBlock;
+        blockReader >> compressedBlock;
+
+        decompressedBlock = qUncompress(compressedBlock);
+        decompressedBlockReader.setBuffer(&decompressedBlock);
+        decompressedBlockReader.open(QIODevice::ReadOnly);
+    }
+
+    return true;
+}
+
+void BlockCompression::close()
+{
+    // flush output buffer
+    if ((openMode() & QIODevice::WriteOnly) && precompressedBlockWriter) {
+        QByteArray compressedBlock = qCompress(precompressedBlockWriter->buffer(), -1);
+        blockWriter << compressedBlock;
+    }
+    basis->close();
+}
+
+void BlockCompression::setBasis(QIODevice *_basis)
+{
+    basis = _basis;
+    blockReader.setDevice(basis);
+    blockWriter.setDevice(basis);
+}
+
+// read from current decompressed block, if out of space, read and decompress another
+// block from basis
+qint64 BlockCompression::readData(char *data, qint64 remaining)
+{
+    qint64 read = 0;
+    while (remaining > 0) {
+        qint64 single_read = decompressedBlockReader.read(data, remaining);
+        if (single_read == -1)
+            qFatal("miss read");
+
+        remaining -= single_read;
+        read += single_read;
+        data += single_read;
+
+        // need a new block
+        if (remaining > 0) {
+            QByteArray compressedBlock;
+            blockReader >> compressedBlock;
+            if (compressedBlock.size() == 0) {
+                return read;
+            }
+            decompressedBlock = qUncompress(compressedBlock);
+
+            decompressedBlockReader.close();
+            decompressedBlockReader.setBuffer(&decompressedBlock);
+            decompressedBlockReader.open(QIODevice::ReadOnly);
+        }
+    }
+    return blockReader.atEnd() && !basis->isReadable() ? -1 : read;
+}
+
+bool BlockCompression::isSequential() const
+{
+    return true;
+}
+
+qint64 BlockCompression::writeData(const char *data, qint64 remaining)
+{
+    qint64 written = 0;
+
+    while (remaining > 0) {
+        // how much more can be put in this buffer?
+        qint64 capacity = blockSize - precompressedBlockWriter->pos();
+
+        // don't try to write beyond capacity 
+        qint64 write_size = qMin(capacity, remaining);
+
+        qint64 singleWrite = precompressedBlockWriter->write(data, write_size);
+        // ignore the error case here, we consdier basis's failure mode the real
+        // end case
+        if (singleWrite == -1)
+            singleWrite = 0;
+
+        remaining -= singleWrite;
+        data += singleWrite;
+        written += singleWrite;
+
+        if (remaining > 0) {
+            QByteArray compressedBlock = qCompress(precompressedBlockWriter->buffer(), -1);
+
+            if (compressedBlock.size() != 0)
+                blockWriter << compressedBlock;
+
+            delete precompressedBlockWriter;
+            precompressedBlockWriter = new QBuffer;
+            precompressedBlockWriter->open(QIODevice::ReadWrite);
+        }
+    }
+    return basis->isWritable() ? written : -1;
+}
+
+
 
 }  // namespace QtUtils
 
diff --git a/openbr/core/qtutils.h b/openbr/core/qtutils.h
index e9cec08..257ccd8 100644
--- a/openbr/core/qtutils.h
+++ b/openbr/core/qtutils.h
@@ -17,6 +17,7 @@
 #ifndef QTUTILS_QTUTILS_H
 #define QTUTILS_QTUTILS_H
 
+#include <QBuffer>
 #include <QByteArray>
 #include <QDir>
 #include <QFile>
@@ -93,6 +94,38 @@ namespace QtUtils
 
     /**** Rect Utilities ****/
     float overlap(const QRectF &r, const QRectF &s);
+
+    
+    class BlockCompression : public QIODevice
+    {
+    public:
+        BlockCompression(QIODevice *_basis);
+        BlockCompression();
+        int blockSize;
+        QIODevice *basis;
+
+        bool open(QIODevice::OpenMode mode);
+
+        void close();
+
+        void setBasis(QIODevice *_basis);
+
+        QDataStream blockReader;
+        QByteArray decompressedBlock;
+        QBuffer decompressedBlockReader;
+
+        // read from current decompressed block, if out of space, read and decompress another
+        // block from basis
+        qint64 readData(char *data, qint64 remaining);
+
+        bool isSequential() const;
+
+        // write to a QByteArray, when max block sized is reached, compress and write
+        // it to basis
+        QBuffer * precompressedBlockWriter;
+        QDataStream blockWriter;
+        qint64 writeData(const char *data, qint64 remaining);
+    };
 }
 
 #endif // QTUTILS_QTUTILS_H
diff --git a/openbr/openbr.cpp b/openbr/openbr.cpp
index b18f019..3162371 100644
--- a/openbr/openbr.cpp
+++ b/openbr/openbr.cpp
@@ -134,9 +134,9 @@ float br_eval_detection(const char *predicted_gallery, const char *truth_gallery
     return EvalDetection(predicted_gallery, truth_gallery, csv, normalize, minSize, maxSize);
 }
 
-float br_eval_landmarking(const char *predicted_gallery, const char *truth_gallery, const char *csv, int normalization_index_a, int normalization_index_b)
+float br_eval_landmarking(const char *predicted_gallery, const char *truth_gallery, const char *csv, int normalization_index_a, int normalization_index_b, int sample_index, int total_examples)
 {
-    return EvalLandmarking(predicted_gallery, truth_gallery, csv, normalization_index_a, normalization_index_b);
+    return EvalLandmarking(predicted_gallery, truth_gallery, csv, normalization_index_a, normalization_index_b, sample_index, total_examples);
 }
 
 void br_eval_regression(const char *predicted_gallery, const char *truth_gallery, const char *predicted_property, const char *truth_property)
diff --git a/openbr/openbr.h b/openbr/openbr.h
index a21c93c..4f32e86 100644
--- a/openbr/openbr.h
+++ b/openbr/openbr.h
@@ -214,8 +214,10 @@ BR_EXPORT float br_eval_detection(const char *predicted_gallery, const char *tru
  * \param csv Optional \c .csv file to contain performance metrics.
  * \param normalization_index_a Optional first index in the list of points to use for normalization.
  * \param normalization_index_b Optional second index in the list of points to use for normalization.
+ * \param sample_index Optional index for sample landmark image in ground truth gallery.
+ * \param total_examples Optional number of accurate and inaccurate examples to display.
  */
-BR_EXPORT float br_eval_landmarking(const char *predicted_gallery, const char *truth_gallery, const char *csv = "", int normalization_index_a = 0, int normalization_index_b = 1);
+BR_EXPORT float br_eval_landmarking(const char *predicted_gallery, const char *truth_gallery, const char *csv = "", int normalization_index_a = 0, int normalization_index_b = 1, int sample_index = 0, int total_examples = 5);
 
 /*!
  * \brief Evaluates regression accuracy to disk.
diff --git a/openbr/plugins/algorithms.cpp b/openbr/plugins/algorithms.cpp
index f04319a..88ecaea 100644
--- a/openbr/plugins/algorithms.cpp
+++ b/openbr/plugins/algorithms.cpp
@@ -31,15 +31,15 @@ class AlgorithmsInitializer : public Initializer
     void initialize() const
     {
         // Face
-        Globals->abbreviations.insert("FaceRecognition", "FaceDetection+Expand+<FaceRecognitionRegistration>+Expand+<FaceRecognitionExtraction>+<FaceRecognitionEmbedding>+<FaceRecognitionQuantization>+SetMetadata(AlgorithmID,-1):MatchProbability(ByteL1)");
-        Globals->abbreviations.insert("GenderClassification", "FaceDetection+Expand+<FaceClassificationRegistration>+Expand+<FaceClassificationExtraction>+<GenderClassifier>+Discard");
-        Globals->abbreviations.insert("AgeRegression", "FaceDetection+Expand+<FaceClassificationRegistration>+Expand+<FaceClassificationExtraction>+<AgeRegressor>+Discard");
+        Globals->abbreviations.insert("FaceRecognition", "FaceDetection+FaceRecognitionRegistration+<FaceRecognitionExtraction>+<FaceRecognitionEmbedding>+<FaceRecognitionQuantization>+SetMetadata(AlgorithmID,-1):Unit(ByteL1)");
+        Globals->abbreviations.insert("GenderClassification", "FaceDetection+Expand+FaceClassificationRegistration+Expand+<FaceClassificationExtraction>+<GenderClassifier>+Discard");
+        Globals->abbreviations.insert("AgeRegression", "FaceDetection+Expand+FaceClassificationRegistration+Expand+<FaceClassificationExtraction>+<AgeRegressor>+Discard");
         Globals->abbreviations.insert("FaceQuality", "Open+Expand+Cascade(FrontalFace)+ASEFEyes+Affine(64,64,0.25,0.35)+ImageQuality+Cvt(Gray)+DFFS+Discard");
         Globals->abbreviations.insert("MedianFace", "Open+Expand+Cascade(FrontalFace)+ASEFEyes+Affine(256,256,0.37,0.45)+Center(Median)");
         Globals->abbreviations.insert("BlurredFaceDetection", "Open+LimitSize(1024)+SkinMask/(Cvt(Gray)+GradientMask)+And+Morph(Erode,16)+LargestConvexArea");
         Globals->abbreviations.insert("DrawFaceDetection", "Open+Cascade(FrontalFace)+Expand+ASEFEyes+Draw(inPlace=true)");
         Globals->abbreviations.insert("ShowFaceDetection", "DrawFaceDetection+Contract+First+Show+Discard");
-        Globals->abbreviations.insert("DownloadFaceRecognition", "Download+Open+ROI+Expand+Cvt(Gray)+Cascade(FrontalFace)+Expand+<FaceRecognitionRegistration>+Expand+<FaceRecognitionExtraction>+<FaceRecognitionEmbedding>+<FaceRecognitionQuantization>+SetMetadata(AlgorithmID,-1):MatchProbability(ByteL1)");
+        Globals->abbreviations.insert("DownloadFaceRecognition", "Download+Open+ROI+Cvt(Gray)+Cascade(FrontalFace)+FaceRecognitionRegistration+<FaceRecognitionExtraction>+<FaceRecognitionEmbedding>+<FaceRecognitionQuantization>+SetMetadata(AlgorithmID,-1):Unit(ByteL1)");
         Globals->abbreviations.insert("OpenBR", "FaceRecognition");
         Globals->abbreviations.insert("GenderEstimation", "GenderClassification");
         Globals->abbreviations.insert("AgeEstimation", "AgeRegression");
@@ -50,7 +50,7 @@ class AlgorithmsInitializer : public Initializer
         // Video
         Globals->abbreviations.insert("DisplayVideo", "FPSLimit(30)+Show(false,[FrameNumber])+Discard");
         Globals->abbreviations.insert("PerFrameDetection", "SaveMat(original)+Cvt(Gray)+Cascade(FrontalFace)+ASEFEyes+RestoreMat(original)+Draw(inPlace=true)+Show(false,[FrameNumber])+Discard");
-        Globals->abbreviations.insert("AgeGenderDemo", "SaveMat(original)+Cvt(Gray)+Cascade(FrontalFace)+Expand+<FaceClassificationRegistration>+<FaceClassificationExtraction>+<AgeRegressor>/<GenderClassifier>+Discard+RestoreMat(original)+Draw(inPlace=true)+DrawPropertiesPoint([Age,Gender],Affine_0,inPlace=true)+SaveMat(original)+Discard+Contract+RestoreMat(original)+FPSCalc+Show(false,[AvgFPS,Age,Gender])+Discard");
+        Globals->abbreviations.insert("AgeGenderDemo", "SaveMat(original)+Cvt(Gray)+Cascade(FrontalFace)+Expand+FaceClassificationRegistration+<FaceClassificationExtraction>+<AgeRegressor>/<GenderClassifier>+Discard+RestoreMat(original)+Draw(inPlace=true)+DrawPropertiesPoint([Age,Gender],Affine_0,inPlace=true)+SaveMat(original)+Discard+Contract+RestoreMat(original)+FPSCalc+Show(false,[AvgFPS,Age,Gender])+Discard");
         Globals->abbreviations.insert("ShowOpticalFlowField", "SaveMat(original)+AggregateFrames(2)+OpticalFlow(useMagnitude=false)+Grid(100,100)+DrawOpticalFlow+FPSLimit(30)+Show(false)+Discard");
         Globals->abbreviations.insert("ShowOpticalFlowMagnitude", "AggregateFrames(2)+OpticalFlow+Normalize(Range,false,0,255)+Cvt(Color)+Draw+FPSLimit(30)+Show(false)+Discard");
         Globals->abbreviations.insert("ShowMotionSegmentation", "DropFrames(5)+AggregateFrames(2)+OpticalFlow+CvtUChar+WatershedSegmentation+DrawSegmentation+Draw+FPSLimit(30)+Show(false)+Discard");
@@ -92,11 +92,11 @@ class AlgorithmsInitializer : public Initializer
         Globals->abbreviations.insert("DenseHOG", "Gradient+RectRegions(8,8,6,6)+Bin(0,360,8)+Hist(8)");
         Globals->abbreviations.insert("DenseSIFT", "(Grid(10,10)+SIFTDescriptor(12)+ByRow)");
         Globals->abbreviations.insert("DenseSIFT2", "(Grid(5,5)+SIFTDescriptor(12)+ByRow)");
-        Globals->abbreviations.insert("FaceRecognitionRegistration", "(ASEFEyes+Affine(88,88,0.25,0.35)+DownsampleTraining(FTE(DFFS),instances=1))");
+        Globals->abbreviations.insert("FaceRecognitionRegistration", "ASEFEyes+Affine(88,88,0.25,0.35)");
         Globals->abbreviations.insert("FaceRecognitionExtraction", "(Mask+DenseSIFT/DenseLBP+DownsampleTraining(PCA(0.95),instances=1)+Normalize(L2)+Cat)");
         Globals->abbreviations.insert("FaceRecognitionEmbedding", "(Dup(12)+RndSubspace(0.05,1)+DownsampleTraining(LDA(0.98),instances=-2)+Cat+DownsampleTraining(PCA(768),instances=1))");
         Globals->abbreviations.insert("FaceRecognitionQuantization", "(Normalize(L1)+Quantize)");
-        Globals->abbreviations.insert("FaceClassificationRegistration", "(ASEFEyes+Affine(56,72,0.33,0.45)+FTE(DFFS))");
+        Globals->abbreviations.insert("FaceClassificationRegistration", "ASEFEyes+Affine(56,72,0.33,0.45)");
         Globals->abbreviations.insert("FaceClassificationExtraction", "((Grid(7,7)+SIFTDescriptor(8)+ByRow)/DenseLBP+DownsampleTraining(PCA(0.95),instances=-1, inputVariable=Gender)+Cat)");
         Globals->abbreviations.insert("AgeRegressor", "DownsampleTraining(Center(Range),instances=-1, inputVariable=Age)+DownsampleTraining(SVM(RBF,EPS_SVR,inputVariable=Age),instances=100, inputVariable=Age)");
         Globals->abbreviations.insert("GenderClassifier", "DownsampleTraining(Center(Range),instances=-1, inputVariable=Gender)+DownsampleTraining(SVM(RBF,C_SVC,inputVariable=Gender),instances=4000, inputVariable=Gender)");
diff --git a/openbr/plugins/cascade.cpp b/openbr/plugins/cascade.cpp
index bd13f7f..a0b282a 100644
--- a/openbr/plugins/cascade.cpp
+++ b/openbr/plugins/cascade.cpp
@@ -252,6 +252,8 @@ class CascadeTransform : public MetaTransform
     void init()
     {
         cascadeResource.setResourceMaker(new CascadeResourceMaker(model));
+        if (model == "Ear" || model == "Eye" || model == "FrontalFace" || model == "ProfileFace")
+            this->trainable = false;
     }
     
     // Train transform
diff --git a/openbr/plugins/draw.cpp b/openbr/plugins/draw.cpp
index 7c6a1a7..3dc6a9b 100644
--- a/openbr/plugins/draw.cpp
+++ b/openbr/plugins/draw.cpp
@@ -42,11 +42,15 @@ class DrawTransform : public UntrainableTransform
     Q_PROPERTY(bool rects READ get_rects WRITE set_rects RESET reset_rects STORED false)
     Q_PROPERTY(bool inPlace READ get_inPlace WRITE set_inPlace RESET reset_inPlace STORED false)
     Q_PROPERTY(int lineThickness READ get_lineThickness WRITE set_lineThickness RESET reset_lineThickness STORED false)
+    Q_PROPERTY(bool named READ get_named WRITE set_named RESET reset_named STORED false)
+    Q_PROPERTY(bool location READ get_location WRITE set_location RESET reset_location STORED false)
     BR_PROPERTY(bool, verbose, false)
     BR_PROPERTY(bool, points, true)
     BR_PROPERTY(bool, rects, true)
     BR_PROPERTY(bool, inPlace, false)
     BR_PROPERTY(int, lineThickness, 1)
+    BR_PROPERTY(bool, named, true)
+    BR_PROPERTY(bool, location, true)
 
     void project(const Template &src, Template &dst) const
     {
@@ -55,11 +59,12 @@ class DrawTransform : public UntrainableTransform
         dst.m() = inPlace ? src.m() : src.m().clone();
 
         if (points) {
-            const QList<Point2f> pointsList = OpenCVUtils::toPoints(src.file.namedPoints() + src.file.points());
+            const QList<Point2f> pointsList = (named) ? OpenCVUtils::toPoints(src.file.points()+src.file.namedPoints()) : OpenCVUtils::toPoints(src.file.points());
             for (int i=0; i<pointsList.size(); i++) {
                 const Point2f &point = pointsList[i];
                 circle(dst, point, 3, color, -1);
-                if (verbose) putText(dst, QString("%1,(%2,%3)").arg(QString::number(i),QString::number(point.x),QString::number(point.y)).toStdString(), point, FONT_HERSHEY_SIMPLEX, 0.5, verboseColor, 1);
+                QString label = (location) ? QString("%1,(%2,%3)").arg(QString::number(i),QString::number(point.x),QString::number(point.y)) : QString("%1").arg(QString::number(i));
+                if (verbose) putText(dst, label.toStdString(), point, FONT_HERSHEY_SIMPLEX, 0.5, verboseColor, 1);
             }
         }
         if (rects) {
diff --git a/openbr/plugins/eigen3.cpp b/openbr/plugins/eigen3.cpp
index 37f3e68..dfe4ba8 100644
--- a/openbr/plugins/eigen3.cpp
+++ b/openbr/plugins/eigen3.cpp
@@ -600,7 +600,7 @@ class SparseLDATransform : public Transform
 
         //Only works on binary class problems for now
         assert(ldaOrig.projection.cols() == 1);
-        float ldaStd = eigStd(ldaOrig.projection);
+        float ldaStd = EigenUtils::stddev(ldaOrig.projection);
         for (int i = 0; i < ldaOrig.projection.rows(); i++)
             if (abs(ldaOrig.projection(i)) > varThreshold * ldaStd)
                 selections.append(i);
diff --git a/openbr/plugins/filter.cpp b/openbr/plugins/filter.cpp
index c57ca1c..48fbe75 100644
--- a/openbr/plugins/filter.cpp
+++ b/openbr/plugins/filter.cpp
@@ -61,11 +61,23 @@ class BlurTransform : public UntrainableTransform
 {
     Q_OBJECT
     Q_PROPERTY(float sigma READ get_sigma WRITE set_sigma RESET reset_sigma STORED false)
+    Q_PROPERTY(bool ROI READ get_ROI WRITE set_ROI RESET reset_ROI STORED false)
     BR_PROPERTY(float, sigma, 1)
+    BR_PROPERTY(bool, ROI, false)
 
     void project(const Template &src, Template &dst) const
     {
-        GaussianBlur(src, dst, Size(0,0), sigma);
+        if (!ROI) GaussianBlur(src, dst, Size(0,0), sigma);
+        else {
+            dst.m() = src.m();
+            foreach (const QRectF &rect, src.file.rects()) {
+                Rect region(rect.x(), rect.y(), rect.width(), rect.height());
+                Mat input = dst.m();
+                Mat output = input.clone();
+                GaussianBlur(input(region), output(region), Size(0,0), sigma);
+                dst.m() = output;
+            }
+        }
     }
 };
 
diff --git a/openbr/plugins/gui.cpp b/openbr/plugins/gui.cpp
index 80f6174..3aab218 100644
--- a/openbr/plugins/gui.cpp
+++ b/openbr/plugins/gui.cpp
@@ -764,7 +764,8 @@ public:
                         }
                         else
                         {
-                            dst[i].file.set(labelSet[idx], rectSet[idx]);
+                            if (labels.isEmpty()) dst[i].file.appendRect(rectSet[idx]);
+                            else dst[i].file.set(labelSet[idx], rectSet[idx]);
                         }
                     }
                 }
diff --git a/openbr/plugins/meta.cpp b/openbr/plugins/meta.cpp
index df880d5..57888a5 100644
--- a/openbr/plugins/meta.cpp
+++ b/openbr/plugins/meta.cpp
@@ -17,6 +17,8 @@
 #include <QFutureSynchronizer>
 #include <QRegularExpression>
 #include <QtConcurrentRun>
+#include <qbuffer.h>
+
 #include "openbr_internal.h"
 #include "openbr/core/common.h"
 #include "openbr/core/opencvutils.h"
@@ -94,17 +96,15 @@ class PipeTransform : public CompositeTransform
 
         int i = 0;
         while (i < transforms.size()) {
-            fprintf(stderr, "\n%s", qPrintable(transforms[i]->objectName()));
-
             // Conditional statement covers likely case that first transform is untrainable
             if (transforms[i]->trainable) {
-                fprintf(stderr, " training...");
+                qDebug() << "Training " << transforms[i]->description() << "\n...";
                 transforms[i]->train(dataLines);
             }
 
             // if the transform is time varying, we can't project it in parallel
             if (transforms[i]->timeVarying()) {
-                fprintf(stderr, "\n%s projecting...", qPrintable(transforms[i]->objectName()));
+                qDebug() << "Projecting " << transforms[i]->description() << "\n...";
                 for (int j=0; j < dataLines.size();j++) {
                     TemplateList junk;
                     splitFTEs(dataLines[j], junk);
@@ -130,7 +130,16 @@ class PipeTransform : public CompositeTransform
                    !transforms[nextTrainableTransform]->timeVarying())
                 nextTrainableTransform++;
 
-            fprintf(stderr, " projecting...");
+            // No more trainable transforms? Don't need any more projects then
+            if (nextTrainableTransform == transforms.size())
+                break;
+
+            fprintf(stderr, "Projecting %s", qPrintable(transforms[i]->description()));
+            for (int j=i+1; j < nextTrainableTransform; j++)
+                fprintf(stderr,"+%s", qPrintable(transforms[j]->description()));
+            fprintf(stderr, "\n...\n");
+            fflush(stderr);
+
             QFutureSynchronizer<void> futures;
             for (int j=0; j < dataLines.size(); j++)
                 futures.addFuture(QtConcurrent::run(this, &PipeTransform::_projectPartial, &dataLines[j], i, nextTrainableTransform));
@@ -510,7 +519,6 @@ class LoadStoreTransform : public MetaTransform
 
 public:
     Transform *transform;
-    QString baseName;
 
     LoadStoreTransform() : transform(NULL) {}
 
@@ -540,8 +548,8 @@ private:
     void init()
     {
         if (transform != NULL) return;
-        if (fileName.isEmpty()) baseName = QRegExp("^[_a-zA-Z0-9]+$").exactMatch(transformString) ? transformString : QtUtils::shortTextHash(transformString);
-        else baseName = fileName;
+        if (fileName.isEmpty()) fileName = QRegExp("^[_a-zA-Z0-9]+$").exactMatch(transformString) ? transformString : QtUtils::shortTextHash(transformString);
+
         if (!tryLoad())
             transform = make(transformString);
         else
@@ -553,19 +561,28 @@ private:
         return transform->timeVarying();
     }
 
-    void train(const TemplateList &data)
+    void train(const QList<TemplateList> &data)
     {
         if (QFileInfo(getFileName()).exists())
             return;
 
         transform->train(data);
 
-        qDebug("Storing %s", qPrintable(baseName));
-        QByteArray byteArray;
-        QDataStream stream(&byteArray, QFile::WriteOnly);
-        stream << transform->description();
+        qDebug("Storing %s", qPrintable(fileName));
+        QtUtils::BlockCompression compressedOut;
+        QFile fout(fileName);
+        QtUtils::touchDir(fout);
+        compressedOut.setBasis(&fout);
+
+        QDataStream stream(&compressedOut);
+        QString desc = transform->description();
+
+        if (!compressedOut.open(QFile::WriteOnly))
+            qFatal("Failed to open %s for writing.", qPrintable(file));
+
+        stream << desc;
         transform->store(stream);
-        QtUtils::writeFile(baseName, byteArray, -1);
+        compressedOut.close();
     }
 
     void project(const Template &src, Template &dst) const
@@ -595,8 +612,8 @@ private:
 
     QString getFileName() const
     {
-        if (QFileInfo(baseName).exists()) return baseName;
-        const QString file = Globals->sdkPath + "/share/openbr/models/transforms/" + baseName;
+        if (QFileInfo(fileName).exists()) return fileName;
+        const QString file = Globals->sdkPath + "/share/openbr/models/transforms/" + fileName;
         return QFileInfo(file).exists() ? file : QString();
     }
 
@@ -606,12 +623,19 @@ private:
         if (file.isEmpty()) return false;
 
         qDebug("Loading %s", qPrintable(file));
-        QByteArray data;
-        QtUtils::readFile(file, data, true);
-        QDataStream stream(&data, QFile::ReadOnly);
+        QFile fin(file);
+        QtUtils::BlockCompression reader(&fin);
+        if (!reader.open(QIODevice::ReadOnly)) {
+            if (QFileInfo(file).exists()) qFatal("Unable to open %s for reading. Check file permissions.", qPrintable(file));
+            else            qFatal("Unable to open %s for reading. File does not exist.", qPrintable(file));
+        }
+
+        QDataStream stream(&reader);
         stream >> transformString;
+
         transform = Transform::make(transformString);
         transform->load(stream);
+
         return true;
     }
 };
diff --git a/openbr/plugins/misc.cpp b/openbr/plugins/misc.cpp
index 71d94b2..af2f793 100644
--- a/openbr/plugins/misc.cpp
+++ b/openbr/plugins/misc.cpp
@@ -115,6 +115,9 @@ private:
     void project(const Template &src, Template &dst) const
     {
         dst.file = src.file;
+        if (Globals->verbose)
+            qDebug("Opening %s", qPrintable(src.file.flat()));
+
         if (src.empty()) {
             const Mat img = imread(src.file.resolved().toStdString(), mode);
             if (img.data) dst.append(img);
@@ -126,6 +129,8 @@ private:
                 else          dst.file.fte = true;
             }
         }
+        if (dst.file.fte)
+            qWarning("Error opening %s", qPrintable(src.file.flat()));
     }
 };
 BR_REGISTER(Transform, ReadTransform)
diff --git a/openbr/plugins/process.cpp b/openbr/plugins/process.cpp
index 548fc11..ba3559b 100644
--- a/openbr/plugins/process.cpp
+++ b/openbr/plugins/process.cpp
@@ -349,6 +349,7 @@ public:
         QByteArray data = readArray;
         QDataStream deserializer(data);
         Transform *res = Transform::deserialize(deserializer);
+        sendSignal(CommunicationManager::OUTPUT_AVAILABLE);
         return res;
     }
 
@@ -534,6 +535,8 @@ struct ProcessData
 class ProcessWrapperTransform : public WrapperTransform
 {
     Q_OBJECT
+    Q_PROPERTY(int concurrentCount READ get_concurrentCount WRITE set_concurrentCount RESET reset_concurrentCount STORED false)
+    BR_PROPERTY(int, concurrentCount, 2)
 
     QString baseKey;
 
@@ -579,17 +582,33 @@ class ProcessWrapperTransform : public WrapperTransform
         if (transform) {
             QDataStream out(&serialized, QFile::WriteOnly);
             transform->serialize(out);
+            tcount = Globals->parallelism;
+            counter.acquire(counter.available());
+            counter.release(this->concurrentCount);
         }
     }
 
-    QByteArray serialized;
+    static QSemaphore counter;
+    mutable int tcount;
+    mutable QByteArray serialized;
     void transmitTForm(CommunicationManager *localComm) const
     {
         if (serialized.isEmpty() )
             qFatal("Trying to transmit empty transform!");
 
+        counter.acquire(1);
+        static QMutex transmission;
+        QMutexLocker lock(&transmission);
+        tcount--;
+
         localComm->writeArray = serialized;
+        if (tcount == 0)
+            serialized.clear();
+        lock.unlock();
+
         emit localComm->pulseSendSerialized();
+        localComm->getSignal();
+        counter.release(1);
     }
 
     void activateProcess(ProcessData *data) const
@@ -633,6 +652,7 @@ public:
     bool processActive;
     ProcessWrapperTransform() : WrapperTransform(false) { processActive = false; }
 };
+QSemaphore ProcessWrapperTransform::counter;
 
 BR_REGISTER(Transform, ProcessWrapperTransform)
 
diff --git a/openbr/plugins/quality.cpp b/openbr/plugins/quality.cpp
index 469009e..11d4bbd 100644
--- a/openbr/plugins/quality.cpp
+++ b/openbr/plugins/quality.cpp
@@ -77,6 +77,12 @@ class ImpostorUniquenessMeasureTransform : public Transform
 
 BR_REGISTER(Transform, ImpostorUniquenessMeasureTransform)
 
+
+float KDEPointer(const QList<float> *scores, double x, double h)
+{
+    return Common::KernelDensityEstimation(*scores, x, h);
+}
+
 /* Kernel Density Estimator */
 struct KDE
 {
@@ -85,20 +91,35 @@ struct KDE
     QList<float> bins;
 
     KDE() : min(0), max(1), mean(0), stddev(1) {}
-    KDE(const QList<float> &scores)
+
+    KDE(const QList<float> &scores, bool trainKDE)
     {
         Common::MinMax(scores, &min, &max);
         Common::MeanStdDev(scores, &mean, &stddev);
+
+        if (!trainKDE)
+            return;
+
         double h = Common::KernelDensityBandwidth(scores);
         const int size = 255;
         bins.reserve(size);
-        for (int i=0; i<size; i++)
-            bins.append(Common::KernelDensityEstimation(scores, min + (max-min)*i/(size-1), h));
+
+        QFutureSynchronizer<float> futures;
+
+        for (int i=0; i < size; i++)
+            futures.addFuture(QtConcurrent::run(KDEPointer, &scores, min + (max-min)*i/(size-1), h));
+        futures.waitForFinished();
+
+        foreach(const QFuture<float> & future, futures.futures())
+            bins.append(future.result());
     }
 
     float operator()(float score, bool gaussian = true) const
     {
         if (gaussian) return 1/(stddev*sqrt(2*CV_PI))*exp(-0.5*pow((score-mean)/stddev, 2));
+        if (bins.empty())
+            return -std::numeric_limits<float>::max();
+	
         if (score <= min) return bins.first();
         if (score >= max) return bins.last();
         const float x = (score-min)/(max-min)*bins.size();
@@ -123,8 +144,8 @@ struct MP
 {
     KDE genuine, impostor;
     MP() {}
-    MP(const QList<float> &genuineScores, const QList<float> &impostorScores)
-        : genuine(genuineScores), impostor(impostorScores) {}
+    MP(const QList<float> &genuineScores, const QList<float> &impostorScores, bool trainKDE)
+        : genuine(genuineScores, trainKDE), impostor(impostorScores, trainKDE) {}
     float operator()(float score, bool gaussian = true) const
     {
         const float g = genuine(score, gaussian);
@@ -165,7 +186,7 @@ class MatchProbabilityDistance : public Distance
         const QList<int> labels = src.indexProperty(inputVariable);
         QScopedPointer<MatrixOutput> matrixOutput(MatrixOutput::make(FileList(src.size()), FileList(src.size())));
         distance->compare(src, src, matrixOutput.data());
-
+	
         QList<float> genuineScores, impostorScores;
         genuineScores.reserve(labels.size());
         impostorScores.reserve(labels.size()*labels.size());
@@ -178,8 +199,8 @@ class MatchProbabilityDistance : public Distance
                 else                        impostorScores.append(score);
             }
         }
-
-        mp = MP(genuineScores, impostorScores);
+	
+        mp = MP(genuineScores, impostorScores, !gaussian);
     }
 
     float compare(const Template &target, const Template &query) const
diff --git a/openbr/plugins/stasm4.cpp b/openbr/plugins/stasm4.cpp
index 0a350e3..98bbad6 100644
--- a/openbr/plugins/stasm4.cpp
+++ b/openbr/plugins/stasm4.cpp
@@ -180,7 +180,7 @@ private:
     void project(const Template &src, Template &dst) const
     {
         QList<float> paramList = src.file.getList<float>("affineParameters");
-        Eigen::MatrixXf points = pointsToMatrix(src.file.points(), true);
+        Eigen::MatrixXf points = EigenUtils::pointsToMatrix(src.file.points(), true);
         Eigen::MatrixXf affine = Eigen::MatrixXf::Zero(3, 3);
         for (int i = 0, cnt = 0; i < 2; i++)
             for (int j = 0; j < 3; j++, cnt++)
@@ -192,7 +192,7 @@ private:
         points =  affineInv * pointsT;
         dst = src;
         dst.file.clearPoints();
-        dst.file.setPoints(matrixToPoints(points.transpose()));
+        dst.file.setPoints(EigenUtils::matrixToPoints(points.transpose()));
     }
 };
 
diff --git a/share/openbr/models b/share/openbr/models
index 79938fe..85842e6 160000
--- a/share/openbr/models
+++ b/share/openbr/models
@@ -1 +1 @@
-Subproject commit 79938fe401faafead086b4711dd0b8f898a7a21e
+Subproject commit 85842e6da7738e317b9d40d5a395b92cd7b996e1