diff --git a/openbr/core/common.h b/openbr/core/common.h
index e95ffe0..9e9b2ae 100644
--- a/openbr/core/common.h
+++ b/openbr/core/common.h
@@ -118,22 +118,22 @@ T Max(const QList<T> &vals)
 /*!
  * \brief Returns the mean and standard deviation of a vector of values.
  */
-template <typename T>
-void Mean(const QList<T> &vals, double *mean)
+template <template<class> class V, typename T>
+void Mean(const V<T> &vals, double *mean)
 {
     const int size = vals.size();
 
     // Compute Mean
     double sum = 0;
-    for (int i=0; i<size; i++) sum += vals[i];
+    foreach (int val, vals) sum += val;
     *mean = (size == 0) ? 0 : sum / size;
 }
 
 /*!
  * \brief Returns the mean and standard deviation of a vector of values.
  */
-template <typename T>
-void MeanStdDev(const QList<T> &vals, double *mean, double *stddev)
+template <template<class> class V, typename T>
+void MeanStdDev(const V<T> &vals, double *mean, double *stddev)
 {
     const int size = vals.size();
 
@@ -141,8 +141,8 @@ void MeanStdDev(const QList<T> &vals, double *mean, double *stddev)
 
     // Compute Standard Deviation
     double variance = 0;
-    for (int i=0; i<size; i++) {
-        double delta = vals[i] - *mean;
+    foreach (T val, vals) {
+        const double delta = val - *mean;
         variance += delta * delta;
     }
     *stddev = (size == 0) ? 0 : sqrt(variance/size);
@@ -193,8 +193,8 @@ QList<T> CumSum(const QList<T> &vals)
 /*!
  * \brief Calculate DKE bandwidth parameter 'h'
  */
-template <typename T>
-double KernelDensityBandwidth(const QList<T> &vals)
+template <template<class> class V, typename T>
+double KernelDensityBandwidth(const V<T> &vals)
 {
     double mean, stddev;
     MeanStdDev(vals, &mean, &stddev);
@@ -204,8 +204,8 @@ double KernelDensityBandwidth(const QList<T> &vals)
 /*!
  * \brief Compute kernel density at value x with bandwidth h.
  */
-template <typename T>
-double KernelDensityEstimation(const QList<T> &vals, double x, double h)
+template <template<class> class V, typename T>
+double KernelDensityEstimation(const V<T> &vals, double x, double h)
 {
     double y = 0;
     foreach (T val, vals)
@@ -320,16 +320,15 @@ QList<T> RemoveOutliers(QList<T> vals)
 /*!
  * \brief Sorts and evenly downsamples a vector to size k.
  */
-template <typename T>
-QList<T> Downsample(QList<T> vals, long k)
+template <template<class> class V, typename T>
+V<T> Downsample(V<T> vals, int k)
 {
-    // Use 'long' instead of 'int' so multiplication doesn't overflow
-    qSort(vals);
-    long size = (long)vals.size();
+    std::sort(vals.begin(), vals.end());
+    int size = vals.size();
     if (size <= k) return vals;
 
-    QList<T> newVals; newVals.reserve(k);
-    for (long i=0; i<k; i++) newVals.push_back(vals[i * (size-1) / (k-1)]);
+    V<T> newVals; newVals.reserve(k);
+    for (int i=0; i<k; i++) newVals.push_back(vals[long(i) * long(size-1) / long(k-1)]);
     return newVals;
 }
 
diff --git a/openbr/plugins/algorithms.cpp b/openbr/plugins/algorithms.cpp
index de277be..2d0f24f 100644
--- a/openbr/plugins/algorithms.cpp
+++ b/openbr/plugins/algorithms.cpp
@@ -42,7 +42,7 @@ class AlgorithmsInitializer : public Initializer
         Globals->abbreviations.insert("OpenBR", "FaceRecognition");
         Globals->abbreviations.insert("GenderEstimation", "GenderClassification");
         Globals->abbreviations.insert("AgeEstimation", "AgeRegression");
-        Globals->abbreviations.insert("FaceRecognitionHoG", "Open+Cvt(Gray)+Cascade(FrontalFace)+ASEFEyes+Affine(64,64,0.25,0.35)+Gradient+Bin(0,360,8,true)+Merge+Integral+IntegralSampler+RootNorm+ProductQuantization(2,true):ProductQuantization(true)");
+        Globals->abbreviations.insert("FaceRecognitionHoG", "Open+Cvt(Gray)+Cascade(FrontalFace)+ASEFEyes+Affine(64,64,0.25,0.35)+Gradient+Bin(0,360,8,true)+Merge+Integral+IntegralSampler+ProductQuantization(2,L1,true):ProductQuantization(true)");
 
         // Generic Image Processing
         Globals->abbreviations.insert("SIFT", "Open+KeyPointDetector(SIFT)+KeyPointDescriptor(SIFT):KeyPointMatcher(BruteForce)");
diff --git a/openbr/plugins/distance.cpp b/openbr/plugins/distance.cpp
index b928427..b987b24 100644
--- a/openbr/plugins/distance.cpp
+++ b/openbr/plugins/distance.cpp
@@ -37,6 +37,7 @@ class DistDistance : public Distance
     Q_OBJECT
     Q_ENUMS(Metric)
     Q_PROPERTY(Metric metric READ get_metric WRITE set_metric RESET reset_metric STORED false)
+    Q_PROPERTY(bool negLogPlusOne READ get_negLogPlusOne WRITE set_negLogPlusOne RESET reset_negLogPlusOne STORED false)
 
 public:
     /*!< */
@@ -51,6 +52,7 @@ public:
 
 private:
     BR_PROPERTY(Metric, metric, L2)
+    BR_PROPERTY(bool, negLogPlusOne, true)
 
     float compare(const Template &a, const Template &b) const
     {
@@ -61,8 +63,7 @@ private:
         float result = std::numeric_limits<float>::max();
         switch (metric) {
           case Correlation:
-            result = -compareHist(a, b, CV_COMP_CORREL);
-            break;
+            return compareHist(a, b, CV_COMP_CORREL);
           case ChiSquared:
             result = compareHist(a, b, CV_COMP_CHISQR);
             break;
@@ -82,8 +83,7 @@ private:
             result = norm(a, b, NORM_L2);
             break;
           case Cosine:
-            result = cosine(a, b);
-            break;
+            return cosine(a, b);
           default:
             qFatal("Invalid metric");
         }
@@ -91,7 +91,7 @@ private:
         if (result != result)
             qFatal("NaN result.");
 
-        return -log(result+1);
+        return negLogPlusOne ? -log(result+1) : result;
     }
 
     static float cosine(const Mat &a, const Mat &b)
diff --git a/openbr/plugins/quantize.cpp b/openbr/plugins/quantize.cpp
index 1302846..8f0b3cc 100644
--- a/openbr/plugins/quantize.cpp
+++ b/openbr/plugins/quantize.cpp
@@ -18,6 +18,7 @@
 #include <QtConcurrentRun>
 #include <openbr/openbr_plugin.h>
 
+#include "openbr/core/common.h"
 #include "openbr/core/opencvutils.h"
 
 using namespace cv;
@@ -133,7 +134,7 @@ class ProductQuantizationDistance : public Distance
             const uchar *bData = b[i].data;
             const float *lut = (const float*)ProductQuantizationLUTs[i].data;
             for (int j=0; j<elements; j++)
-                 distance += lut[i*256*256 + aData[j]*256+bData[j]];
+                 distance += lut[j*256*256 + aData[j]*256+bData[j]];
         }
         if (!bayesian) distance = -log(distance+1);
         return distance;
@@ -151,8 +152,10 @@ class ProductQuantizationTransform : public Transform
 {
     Q_OBJECT
     Q_PROPERTY(int n READ get_n WRITE set_n RESET reset_n STORED false)
+    Q_PROPERTY(br::Distance *distance READ get_distance WRITE set_distance RESET reset_distance STORED false)
     Q_PROPERTY(bool bayesian READ get_bayesian WRITE set_bayesian RESET reset_bayesian STORED false)
     BR_PROPERTY(int, n, 2)
+    BR_PROPERTY(br::Distance*, distance, Distance::make("L2", this))
     BR_PROPERTY(bool, bayesian, false)
 
     int index;
@@ -166,76 +169,82 @@ public:
     }
 
 private:
-    static double likelihoodRatio(const QPair<int,int> &totals, const QList<int> &targets, const QList<int> &queries)
+    void _train(const Mat &data, const QList<int> &labels, Mat *lut, Mat *center)
     {
-        int positives = 1, negatives = 1; // Equal priors
-        foreach (int t, targets)
-            foreach (int q, queries)
-                if (t == q) positives++;
-                else        negatives++;
-        return log((float(positives)/float(totals.first)) / (float(negatives)/float(totals.second)));
-    }
+        Mat clusterLabels;
+        kmeans(data, 256, clusterLabels, TermCriteria(TermCriteria::MAX_ITER, 10, 0), 3, KMEANS_PP_CENTERS, *center);
 
-    void _train(const Mat &data, const QPair<int,int> &totals, Mat &lut, int i, const QList<int> &templateLabels)
-    {
-        Mat labels, center;
-        kmeans(data.colRange(i*n,(i+1)*n), 256, labels, TermCriteria(TermCriteria::MAX_ITER, 10, 0), 3, KMEANS_PP_CENTERS, center);
-        QList<int> clusterLabels = OpenCVUtils::matrixToVector<int>(labels);
+        for (int j=0; j<256; j++)
+            for (int k=0; k<256; k++)
+                lut->at<float>(0,j*256+k) = distance->compare(center->row(j), center->row(k));
+
+        if (!bayesian) return;
+
+        QList<int> indicies = OpenCVUtils::matrixToVector<int>(clusterLabels);
+        QVector<float> genuineScores; genuineScores.reserve(data.rows);
+        QVector<float> impostorScores; impostorScores.reserve(data.rows*data.rows/2);
+        for (int i=0; i<indicies.size(); i++)
+            for (int j=i+1; j<indicies.size(); j++) {
+                const float score = lut->at<float>(0, indicies[i]*256+indicies[j]);
+                if (labels[i] == labels[j]) genuineScores.append(score);
+                else                        impostorScores.append(score);
+            }
+        genuineScores = Common::Downsample(genuineScores, 256);
+        impostorScores = Common::Downsample(impostorScores, 256);
 
-        QHash< int, QList<int> > clusters; // QHash<clusterLabel, QList<templateLabel>>
-        for (int j=0; j<clusterLabels.size(); j++)
-            clusters[clusterLabels[j]].append(templateLabels[j]);
+        double hGenuine = Common::KernelDensityBandwidth(genuineScores);
+        double hImpostor = Common::KernelDensityBandwidth(impostorScores);
 
         for (int j=0; j<256; j++)
             for (int k=0; k<256; k++)
-                lut.at<float>(i,j*256+k) = bayesian ? likelihoodRatio(totals, clusters[j], clusters[k]) :
-                                                      norm(center.row(j), center.row(k), NORM_L2);
-        centers[i] = center;
+                lut->at<float>(0,j*256+k) = log(Common::KernelDensityEstimation(genuineScores, lut->at<float>(0,j*256+k), hGenuine) /
+                                                Common::KernelDensityEstimation(impostorScores, lut->at<float>(0,j*256+k), hImpostor));
     }
 
     void train(const TemplateList &src)
     {
         Mat data = OpenCVUtils::toMat(src.data());
         if (data.cols % n != 0) qFatal("Expected dimensionality to be divisible by n.");
-        const QList<int> templateLabels = src.labels<int>();
-        int totalPositives = 0, totalNegatives = 0;
-        for (int i=0; i<templateLabels.size(); i++)
-            for (int j=0; j<templateLabels.size(); j++)
-                if (templateLabels[i] == templateLabels[j]) totalPositives++;
-                else                                        totalNegatives++;
-        QPair<int,int> totals(totalPositives, totalNegatives);
+        const QList<int> labels = src.labels<int>();
 
         Mat &lut = ProductQuantizationLUTs[index];
         lut = Mat(data.cols/n, 256*256, CV_32FC1);
 
-        for (int i=0; i<lut.rows; i++)
+        QList<Mat> subdata, subluts;
+        for (int i=0; i<lut.rows; i++) {
             centers.append(Mat());
+            subdata.append(data.colRange(i*n,(i+1)*n));
+            subluts.append(lut.row(i));
+        }
 
         QFutureSynchronizer<void> futures;
         for (int i=0; i<lut.rows; i++) {
-            if (Globals->parallelism) futures.addFuture(QtConcurrent::run(this, &ProductQuantizationTransform::_train, data, totals, lut, i, templateLabels));
-            else                                                                                               _train (data, totals, lut, i, templateLabels);
+            if (Globals->parallelism) futures.addFuture(QtConcurrent::run(this, &ProductQuantizationTransform::_train, subdata[i], labels, &subluts[i], &centers[i]));
+            else                                                                                               _train (subdata[i], labels, &subluts[i], &centers[i]);
         }
         futures.waitForFinished();
     }
 
+    int getIndex(const Mat &m, const Mat &center) const
+    {
+        int bestIndex = 0;
+        double bestDistance = std::numeric_limits<double>::max();
+        for (int j=0; j<256; j++) {
+            double distance = norm(m, center.row(j), NORM_L2);
+            if (distance < bestDistance) {
+                bestDistance = distance;
+                bestIndex = j;
+            }
+        }
+        return bestIndex;
+    }
+
     void project(const Template &src, Template &dst) const
     {
         Mat m = src.m().reshape(1, 1);
         dst = Mat(1, m.cols/n, CV_8UC1);
-        for (int i=0; i<dst.m().cols; i++) {
-            int bestIndex = 0;
-            double bestDistance = std::numeric_limits<double>::max();
-            Mat m_i = m.colRange(i*n, (i+1)*n);
-            for (int j=0; j<256; j++) {
-                double distance = norm(m_i, centers[index].row(j), NORM_L2);
-                if (distance < bestDistance) {
-                    bestDistance = distance;
-                    bestIndex = j;
-                }
-            }
-            dst.m().at<uchar>(0,i) = bestIndex;
-        }
+        for (int i=0; i<dst.m().cols; i++)
+            dst.m().at<uchar>(0,i) = getIndex(m.colRange(i*n, (i+1)*n), centers[i]);
     }
 
     void store(QDataStream &stream) const