diff --git a/openbr/core/bee.cpp b/openbr/core/bee.cpp
index 05143df..e0ede3f 100644
--- a/openbr/core/bee.cpp
+++ b/openbr/core/bee.cpp
@@ -234,22 +234,26 @@ void BEE::writeMask(const Mat &m, const QString &mask, const QString &targetSigs
 void BEE::makeMask(const QString &targetInput, const QString &queryInput, const QString &mask)
 {
     qDebug("Making mask from %s and %s to %s", qPrintable(targetInput), qPrintable(queryInput), qPrintable(mask));
+    FileList targes = TemplateList::fromGallery(targetInput).files();
+    FileList queries = (queryInput == ".") ? targes : TemplateList::fromGallery(queryInput).files();
+    writeMask(makeMask(targes, queries), mask, targetInput, queryInput);
+}
 
-    FileList targetFiles = TemplateList::fromGallery(targetInput).files();
-    FileList queryFiles = (queryInput == ".") ? targetFiles : TemplateList::fromGallery(queryInput).files();
-    QList<float> targetLabels = targetFiles.labels();
-    QList<float> queryLabels = queryFiles.labels();
-    QList<int> targetPartitions = targetFiles.crossValidationPartitions();
-    QList<int> queryPartitions = queryFiles.crossValidationPartitions();
-
-    Mat vals(queryFiles.size(), targetFiles.size(), CV_8UC1);
-    for (int i=0; i<queryFiles.size(); i++) {
-        const QString &fileA = queryFiles[i];
+cv::Mat BEE::makeMask(const br::FileList &targets, const br::FileList &queries)
+{
+    QList<float> targetLabels = targets.labels();
+    QList<float> queryLabels = queries.labels();
+    QList<int> targetPartitions = targets.crossValidationPartitions();
+    QList<int> queryPartitions = queries.crossValidationPartitions();
+
+    Mat mask(queries.size(), targets.size(), CV_8UC1);
+    for (int i=0; i<queries.size(); i++) {
+        const QString &fileA = queries[i];
         const int labelA = queryLabels[i];
         const int partitionA = queryPartitions[i];
 
-        for (int j=0; j<targetFiles.size(); j++) {
-            const QString &fileB = targetFiles[j];
+        for (int j=0; j<targets.size(); j++) {
+            const QString &fileB = targets[j];
             const int labelB = targetLabels[j];
             const int partitionB = targetPartitions[j];
 
@@ -260,10 +264,11 @@ void BEE::makeMask(const QString &targetInput, const QString &queryInput, const 
             else if (partitionA != partitionB) val = DontCare;
             else if (labelA == labelB)         val = Match;
             else                               val = NonMatch;
-            vals.at<Mask_t>(i,j) = val;
+            mask.at<Mask_t>(i,j) = val;
         }
     }
-    writeMask(vals, mask, targetInput, queryInput);
+
+    return mask;
 }
 
 void BEE::combineMasks(const QStringList &inputMasks, const QString &outputMask, const QString &method)
diff --git a/openbr/core/bee.h b/openbr/core/bee.h
index 6880944..3cd9fcd 100644
--- a/openbr/core/bee.h
+++ b/openbr/core/bee.h
@@ -46,6 +46,7 @@ namespace BEE
 
     // Write BEE files
     void makeMask(const QString &targetInput, const QString &queryInput, const QString &mask);
+    cv::Mat makeMask(const br::FileList &targets, const br::FileList &queries);
     void combineMasks(const QStringList &inputMasks, const QString &outputMask, const QString &method);
 }
 
diff --git a/openbr/core/common.h b/openbr/core/common.h
index 9e9b2ae..4bb9620 100644
--- a/openbr/core/common.h
+++ b/openbr/core/common.h
@@ -119,14 +119,14 @@ T Max(const QList<T> &vals)
  * \brief Returns the mean and standard deviation of a vector of values.
  */
 template <template<class> class V, typename T>
-void Mean(const V<T> &vals, double *mean)
+double Mean(const V<T> &vals)
 {
     const int size = vals.size();
 
     // Compute Mean
     double sum = 0;
     foreach (int val, vals) sum += val;
-    *mean = (size == 0) ? 0 : sum / size;
+    return (size == 0) ? 0 : sum / size;
 }
 
 /*!
@@ -137,7 +137,7 @@ void MeanStdDev(const V<T> &vals, double *mean, double *stddev)
 {
     const int size = vals.size();
 
-    Mean(vals, mean);
+    *mean = Mean(vals);
 
     // Compute Standard Deviation
     double variance = 0;
diff --git a/openbr/core/plot.cpp b/openbr/core/plot.cpp
index c3056fb..03b8005 100644
--- a/openbr/core/plot.cpp
+++ b/openbr/core/plot.cpp
@@ -113,9 +113,6 @@ float Evaluate(const QString &simmat, const QString &mask, const QString &csv)
 {
     qDebug("Evaluating %s with %s", qPrintable(simmat), qPrintable(mask));
 
-    const int Max_Points = 500;
-    float result = -1;
-
     // Read files
     const Mat scores = BEE::readSimmat(simmat);
     File maskFile(mask);
@@ -124,13 +121,21 @@ float Evaluate(const QString &simmat, const QString &mask, const QString &csv)
     const Mat masks = BEE::readMask(maskFile);
     if (scores.size() != masks.size()) qFatal("Simmat (%i,%i) / Mask (%i,%i) size mismatch.", scores.rows, scores.cols, masks.rows, masks.cols);
 
+    return Evaluate(scores, masks, csv);
+}
+
+float Evaluate(const Mat &simmat, const Mat &mask, const QString &csv)
+{
+    const int Max_Points = 500;
+    float result = -1;
+
     // Make comparisons
-    QList<Comparison> comparisons; comparisons.reserve(scores.rows*scores.cols);
+    QList<Comparison> comparisons; comparisons.reserve(simmat.rows*simmat.cols);
     int genuineCount = 0, impostorCount = 0, numNaNs = 0;
-    for (int i=0; i<scores.rows; i++) {
-        for (int j=0; j<scores.cols; j++) {
-            const BEE::Mask_t mask_val = masks.at<BEE::Mask_t>(i,j);
-            const BEE::Simmat_t simmat_val = scores.at<BEE::Simmat_t>(i,j);
+    for (int i=0; i<simmat.rows; i++) {
+        for (int j=0; j<simmat.cols; j++) {
+            const BEE::Mask_t mask_val = mask.at<BEE::Mask_t>(i,j);
+            const BEE::Simmat_t simmat_val = simmat.at<BEE::Simmat_t>(i,j);
             if (mask_val == BEE::DontCare) continue;
             if (simmat_val != simmat_val) { numNaNs++; continue; }
             comparisons.append(Comparison(simmat_val, j, i, mask_val == BEE::Match));
@@ -149,7 +154,7 @@ float Evaluate(const QString &simmat, const QString &mask, const QString &csv)
     double genuineSum = 0, impostorSum = 0;
     QList<OperatingPoint> operatingPoints;
     QList<float> genuines, impostors;
-    QVector<int> firstGenuineReturns(scores.rows, 0);
+    QVector<int> firstGenuineReturns(simmat.rows, 0);
 
     int falsePositives = 0, previousFalsePositives = 0;
     int truePositives = 0, previousTruePositives = 0;
@@ -202,11 +207,11 @@ float Evaluate(const QString &simmat, const QString &mask, const QString &csv)
     // Write Metadata table
     QStringList lines;
     lines.append("Plot,X,Y");
-    lines.append("Metadata,"+QString::number(scores.cols)+",Gallery");
-    lines.append("Metadata,"+QString::number(scores.rows)+",Probe");
+    lines.append("Metadata,"+QString::number(simmat.cols)+",Gallery");
+    lines.append("Metadata,"+QString::number(simmat.rows)+",Probe");
     lines.append("Metadata,"+QString::number(genuineCount)+",Genuine");
     lines.append("Metadata,"+QString::number(impostorCount)+",Impostor");
-    lines.append("Metadata,"+QString::number(scores.cols*scores.rows-(genuineCount+impostorCount))+",Ignored");
+    lines.append("Metadata,"+QString::number(simmat.cols*simmat.rows-(genuineCount+impostorCount))+",Ignored");
 
     // Write Detection Error Tradeoff (DET), PRE, REC
     int points = qMin(operatingPoints.size(), Max_Points);
@@ -382,7 +387,11 @@ struct RPlot
                        "ERR$Y <- as.numeric(as.character(ERR$Y))\n"
                        "SD$Y <- as.factor(unique(as.character(SD$Y)))\n"
                        "BC$Y <- as.numeric(as.character(BC$Y))\n"
-                       "CMC$Y <- as.numeric(as.character(CMC$Y))\n");
+                       "CMC$Y <- as.numeric(as.character(CMC$Y))\n"
+                       "\n"
+                       "# Code to format FAR values\n"
+                       "far_names <- list('0.001'=\"FAR = 0.1%\", '0.01'=\"FAR = 1%\")\n"
+                       "far_labeller <- function(variable,value) { return(far_names[as.character(value)]) }\n");
 
         // Open output device
         file.write(qPrintable(QString("\n"
@@ -461,7 +470,7 @@ bool Plot(const QStringList &files, const br::File &destination, bool show)
                             QString(", xlab=\"False Accept Rate\", ylab=\"True Accept Rate\") + theme_minimal()") +
                             (p.major.size > 1 ? getScale("colour", p.major.header, p.major.size) : QString()) +
                             (p.minor.size > 1 ? QString(" + scale_linetype_discrete(\"%1\")").arg(p.minor.header) : QString()) +
-                            QString(" + scale_x_log10() + scale_y_continuous(labels=percent)") +
+                            QString(" + scale_x_log10(labels=percent) + scale_y_continuous(labels=percent) + annotation_logticks(sides=\"b\")") +
                             QString("\nggsave(\"%1\")\n").arg(p.subfile("ROC"))));
 
     p.file.write(qPrintable(QString("qplot(X, Y, data=DET%1").arg((p.major.smooth || p.minor.smooth) ? ", geom=\"smooth\", method=loess, level=0.99" : ", geom=\"line\"") +
@@ -470,7 +479,7 @@ bool Plot(const QStringList &files, const br::File &destination, bool show)
                             QString(", xlab=\"False Accept Rate\", ylab=\"False Reject Rate\") + geom_abline(alpha=0.5, colour=\"grey\", linetype=\"dashed\") + theme_minimal()") +
                             (p.major.size > 1 ? getScale("colour", p.major.header, p.major.size) : QString()) +
                             (p.minor.size > 1 ? QString(" + scale_linetype_discrete(\"%1\")").arg(p.minor.header) : QString()) +
-                            QString(" + scale_x_log10() + scale_y_log10()") +
+                            QString(" + scale_x_log10(labels=percent) + scale_y_log10(labels=percent) + annotation_logticks()") +
                             QString("\nggsave(\"%1\")\n").arg(p.subfile("DET"))));
 
     p.file.write(qPrintable(QString("qplot(X, data=SD, geom=\"histogram\", fill=Y, position=\"identity\", alpha=I(1/2)") +
@@ -495,7 +504,7 @@ bool Plot(const QStringList &files, const br::File &destination, bool show)
                             QString(", xlab=\"%1False Accept Rate\"").arg(p.major.size > 1 ? p.major.header + " / " : QString()) +
                             QString(", ylab=\"True Accept Rate%1\") + theme_minimal()").arg(p.minor.size > 1 ? " / " + p.minor.header : QString()) +
                             (p.major.size > 1 ? getScale("fill", p.major.header, p.major.size) : QString()) +
-                            (p.minor.size > 1 ? QString(" + facet_grid(%2 ~ X)").arg(p.minor.header) : QString(" + facet_wrap(~ X)")) +
+                            (p.minor.size > 1 ? QString(" + facet_grid(%2 ~ X)").arg(p.minor.header) : QString(" + facet_grid(. ~ X, labeller=far_labeller)")) +
                             QString(" + scale_y_continuous(labels=percent) + theme(legend.position=\"none\", axis.text.x=element_text(angle=-90, hjust=0))%1").arg((p.major.smooth || p.minor.smooth) ? "" : " + geom_text(data=BC, aes(label=Y, y=0.05))") +
                             QString("\nggsave(\"%1\")\n").arg(p.subfile("BC"))));
 
@@ -503,7 +512,7 @@ bool Plot(const QStringList &files, const br::File &destination, bool show)
                             ((p.flip ? p.major.size : p.minor.size) > 1 ? QString(", colour=factor(%1)").arg(p.flip ? p.major.header : p.minor.header) : QString()) +
                             QString(", xlab=\"Score%1\", ylab=\"Error Rate\") + theme_minimal()").arg((p.flip ? p.minor.size : p.major.size) > 1 ? " / " + (p.flip ? p.minor.header : p.major.header) : QString()) +
                             ((p.flip ? p.major.size : p.minor.size) > 1 ? getScale("colour", p.flip ? p.major.header : p.minor.header, p.flip ? p.major.size : p.minor.size) : QString()) +
-                            QString(" + scale_y_log10()") +
+                            QString(" + scale_y_log10(labels=percent) + annotation_logticks(sides=\"l\")") +
                             ((p.flip ? p.minor.size : p.major.size) > 1 ? QString(" + facet_wrap(~ %1, scales=\"free_x\")").arg(p.flip ? p.minor.header : p.major.header) : QString()) +
                             QString(" + theme(aspect.ratio=1)") +
                             QString("\nggsave(\"%1\")\n").arg(p.subfile("ERR"))));
diff --git a/openbr/core/plot.h b/openbr/core/plot.h
index 782de18..20486e7 100644
--- a/openbr/core/plot.h
+++ b/openbr/core/plot.h
@@ -27,6 +27,7 @@ namespace br
 
 void Confusion(const QString &file, float score, int &true_positives, int &false_positives, int &true_negatives, int &false_negatives);
 float Evaluate(const QString &simmat, const QString &mask, const QString &csv = ""); // Returns TAR @ FAR = 0.01
+float Evaluate(const cv::Mat &scores, const cv::Mat &masks, const QString &csv = "");
 bool Plot(const QStringList &files, const br::File &destination, bool show = false);
 bool PlotMetadata(const QStringList &files, const QString &destination, bool show = false);
 
diff --git a/openbr/plugins/crop.cpp b/openbr/plugins/crop.cpp
index 456cc50..e9cfff4 100644
--- a/openbr/plugins/crop.cpp
+++ b/openbr/plugins/crop.cpp
@@ -90,6 +90,25 @@ BR_REGISTER(Transform, LimitSizeTransform)
 
 /*!
  * \ingroup transforms
+ * \brief Enforce a multiple of \em n columns.
+ * \author Josh Klontz \cite jklontz
+ */
+class DivTransform : public UntrainableTransform
+{
+    Q_OBJECT
+    Q_PROPERTY(int n READ get_n WRITE set_n RESET reset_n STORED false)
+    BR_PROPERTY(int, n, 1)
+
+    void project(const Template &src, Template &dst) const
+    {
+        dst = Mat(src, Rect(0,0,n*(src.m().cols/n),src.m().rows));
+    }
+};
+
+BR_REGISTER(Transform, DivTransform)
+
+/*!
+ * \ingroup transforms
  * \brief Crop out black borders
  * \author Josh Klontz \cite jklontz
  */
diff --git a/openbr/plugins/filter.cpp b/openbr/plugins/filter.cpp
index 94347cb..32e9ec5 100644
--- a/openbr/plugins/filter.cpp
+++ b/openbr/plugins/filter.cpp
@@ -235,7 +235,7 @@ class PowTransform : public UntrainableTransform
 
     void project(const Template &src, Template &dst) const
     {
-        pow(src, power, dst);
+        pow(preserveSign ? abs(src) : src.m(), power, dst);
         if (preserveSign) subtract(Scalar::all(0), dst, dst, src.m() < 0);
     }
 };
diff --git a/openbr/plugins/integral.cpp b/openbr/plugins/integral.cpp
index f6978d4..43fe4f8 100644
--- a/openbr/plugins/integral.cpp
+++ b/openbr/plugins/integral.cpp
@@ -38,15 +38,19 @@ class IntegralSamplerTransform : public UntrainableTransform
     Q_PROPERTY(float scaleFactor READ get_scaleFactor WRITE set_scaleFactor RESET reset_scaleFactor STORED false)
     Q_PROPERTY(float stepFactor READ get_stepFactor WRITE set_stepFactor RESET reset_stepFactor STORED false)
     Q_PROPERTY(int minSize READ get_minSize WRITE set_minSize RESET reset_minSize STORED false)
+    Q_PROPERTY(bool secondOrder READ get_secondOrder WRITE set_secondOrder RESET reset_secondOrder STORED false)
     BR_PROPERTY(int, scales, 6)
     BR_PROPERTY(float, scaleFactor, 2)
     BR_PROPERTY(float, stepFactor, 0.75)
-    BR_PROPERTY(int, minSize, 6)
+    BR_PROPERTY(int, minSize, 8)
+    BR_PROPERTY(bool, secondOrder, false)
 
     void project(const Template &src, Template &dst) const
     {
         typedef Eigen::Map< const Eigen::Matrix<qint32,Eigen::Dynamic,1> > InputDescriptor;
+        typedef Eigen::Map< const Eigen::Matrix<float,Eigen::Dynamic,1> > SecondOrderInputDescriptor;
         typedef Eigen::Map< Eigen::Matrix<float,Eigen::Dynamic,1> > OutputDescriptor;
+
         const Mat &m = src.m();
         if (m.depth() != CV_32S) qFatal("Expected CV_32S matrix depth.");
         const int channels = m.channels();
@@ -56,8 +60,10 @@ class IntegralSamplerTransform : public UntrainableTransform
         float idealSize = min(m.rows, m.cols)-1;
         for (int scale=0; scale<scales; scale++) {
             const int currentSize(idealSize);
-            descriptors += (1+(m.rows-currentSize-1)/int(idealSize*stepFactor)) *
-                           (1+(m.cols-currentSize-1)/int(idealSize*stepFactor));
+            const int numDown = 1+(m.rows-currentSize-1)/int(idealSize*stepFactor);
+            const int numAcross = 1+(m.cols-currentSize-1)/int(idealSize*stepFactor);
+            descriptors += numDown*numAcross;
+            if (secondOrder) descriptors += numDown*(numAcross-1) + (numDown-1)*numAcross;
             idealSize /= scaleFactor;
             if (idealSize < minSize) break;
         }
@@ -81,6 +87,26 @@ class IntegralSamplerTransform : public UntrainableTransform
                     index++;
                 }
             }
+            if (secondOrder) {
+                const int numDown = 1+(m.rows-currentSize-1)/currentStep;
+                const int numAcross = 1+(m.cols-currentSize-1)/currentStep;
+                const float *dataIn = n.ptr<float>(index - numDown*numAcross);
+                for (int i=0; i<numDown; i++) {
+                    for (int j=0; j<numAcross; j++) {
+                        SecondOrderInputDescriptor a(dataIn + (i*numAcross+j)*channels, channels, 1);
+                        if (j < numAcross-1) {
+                            OutputDescriptor y(dataOut+(index*channels), channels, 1);
+                            y = a - SecondOrderInputDescriptor(dataIn + (i*numAcross+j+1)*channels, channels, 1);
+                            index++;
+                        }
+                        if (i < numDown-1) {
+                            OutputDescriptor y(dataOut+(index*channels), channels, 1);
+                            y = a - SecondOrderInputDescriptor(dataIn + ((i+1)*numAcross+j)*channels, channels, 1);
+                            index++;
+                        }
+                    }
+                }
+            }
             idealSize /= scaleFactor;
             if (idealSize < minSize) break;
         }
diff --git a/openbr/plugins/meta.cpp b/openbr/plugins/meta.cpp
index f19e4cc..4c7eee1 100644
--- a/openbr/plugins/meta.cpp
+++ b/openbr/plugins/meta.cpp
@@ -418,6 +418,8 @@ public:
 private:
     void init()
     {
+        if (!transform) return;
+
         trainable = transform->trainable;
         if (!cache.isEmpty()) return;
 
diff --git a/openbr/plugins/output.cpp b/openbr/plugins/output.cpp
index ecd9250..d431fd1 100644
--- a/openbr/plugins/output.cpp
+++ b/openbr/plugins/output.cpp
@@ -40,6 +40,7 @@
 #include "openbr/core/bee.h"
 #include "openbr/core/common.h"
 #include "openbr/core/opencvutils.h"
+#include "openbr/core/plot.h"
 #include "openbr/core/qtutils.h"
 
 namespace br
@@ -259,6 +260,23 @@ BR_REGISTER(Output, EmptyOutput)
 
 /*!
  * \ingroup outputs
+ * \brief Evaluate the output matrix.
+ * \author Josh Klontz \cite jklontz
+ */
+class evalOutput : public MatrixOutput
+{
+    Q_OBJECT
+
+    ~evalOutput()
+    {
+        Evaluate(data, BEE::makeMask(targetFiles, queryFiles), "");
+    }
+};
+
+BR_REGISTER(Output, evalOutput)
+
+/*!
+ * \ingroup outputs
  * \brief Outputs highest ranked matches with scores.
  * \author Scott Klum \cite sklum
  */
diff --git a/openbr/plugins/quality.cpp b/openbr/plugins/quality.cpp
index 3187e95..5ea1c3e 100644
--- a/openbr/plugins/quality.cpp
+++ b/openbr/plugins/quality.cpp
@@ -1,6 +1,9 @@
+#include <QFutureSynchronizer>
+#include <QtConcurrent>
 #include <openbr/openbr_plugin.h>
 
 #include "openbr/core/common.h"
+#include "openbr/core/opencvutils.h"
 
 namespace br
 {
@@ -32,8 +35,7 @@ class ImpostorUniquenessMeasureTransform : public Transform
         QList<float> scores = distance->compare(subset, probe);
         float min, max;
         Common::MinMax(scores, &min, &max);
-        double mean;
-        Common::Mean(scores, &mean);
+        double mean = Common::Mean(scores);
         return (max-mean)/(max-min);
     }
 
@@ -199,6 +201,89 @@ class MatchProbabilityDistance : public Distance
 BR_REGISTER(Distance, MatchProbabilityDistance)
 
 /*!
+ * \ingroup transforms
+ * \brief Normalize by Bhattacharyya coefficient.
+ * \author Josh Klontz \cite jklontz
+ */
+class BhattacharyyaTransform : public Transform
+{
+    Q_OBJECT
+    cv::Mat scale;
+
+    static float bhattacharyyaCoefficient(const cv::Mat &data, const QList<int> &labels)
+    {
+        const QList<float> vals = OpenCVUtils::matrixToVector<float>(data);
+        if (vals.size() != labels.size())
+            qFatal("Logic error.");
+
+        QList<float> genuineScores; genuineScores.reserve(vals.size());
+        QList<float> impostorScores; impostorScores.reserve(vals.size()*vals.size()/2);
+        for (int i=0; i<vals.size(); i++)
+            for (int j=i+1; j<vals.size(); j++)
+                if (labels[i] == labels[j]) genuineScores.append(fabs(vals[i]-vals[j]));
+                else                        impostorScores.append(fabs(vals[i]-vals[j]));
+
+        genuineScores = Common::Downsample(genuineScores, 256);
+        impostorScores = Common::Downsample(impostorScores, 256);
+        double hGenuine = Common::KernelDensityBandwidth(genuineScores);
+        double hImpostor = Common::KernelDensityBandwidth(impostorScores);
+
+        float genuineMin, genuineMax, impostorMin, impostorMax, min, max;
+        Common::MinMax(genuineScores, &genuineMin, &genuineMax);
+        Common::MinMax(impostorScores, &impostorMin, &impostorMax);
+        min = std::min(genuineMin, impostorMin);
+        max = std::max(genuineMax, impostorMax);
+
+        const int steps = 512;
+        double bc = 0;
+        for (int i=0; i<steps; i++) {
+            float score = min + i*(max-min)/(steps-1);
+            bc += sqrt(Common::KernelDensityEstimation(genuineScores, score, hGenuine) *
+                       Common::KernelDensityEstimation(impostorScores, score, hImpostor))/steps;
+        }
+        if ((bc <= 0) || (bc > 1)) qFatal("Logic error.");
+        return -log(bc);
+    }
+
+    void train(const TemplateList &src)
+    {
+        const cv::Mat data = OpenCVUtils::toMat(src.data());
+        const QList<int> labels = src.labels<int>();
+
+        QFutureSynchronizer<float> futures;
+        QList<float> coefficients; coefficients.reserve(data.cols);
+        for (int i=0; i<data.cols; i++)
+            if (Globals->parallelism) futures.addFuture(QtConcurrent::run(&BhattacharyyaTransform::bhattacharyyaCoefficient, data.col(i), labels));
+            else                      coefficients.append(                                         bhattacharyyaCoefficient( data.col(i), labels));
+        futures.waitForFinished();
+        foreach (const QFuture<float> &future, futures.futures())
+            coefficients.append(future.result());
+
+        scale = cv::Mat(1, coefficients.size(), CV_32FC1);
+        for (int i=0; i<coefficients.size(); i++)
+            scale.at<float>(0,i) = coefficients[i];
+
+    }
+
+    void project(const Template &src, Template &dst) const
+    {
+        cv::multiply(src.m().reshape(1,1), scale, dst);
+    }
+
+    void store(QDataStream &stream) const
+    {
+        stream << scale;
+    }
+
+    void load(QDataStream &stream)
+    {
+        stream >> scale;
+    }
+};
+
+BR_REGISTER(Transform, BhattacharyyaTransform)
+
+/*!
  * \ingroup distances
  * \brief Linear normalizes of a distance so the mean impostor score is 0 and the mean genuine score is 1.
  * \author Josh Klontz \cite jklontz
diff --git a/openbr/plugins/quantize.cpp b/openbr/plugins/quantize.cpp
index 8f0b3cc..a6661d3 100644
--- a/openbr/plugins/quantize.cpp
+++ b/openbr/plugins/quantize.cpp
@@ -56,6 +56,72 @@ class QuantizeTransform : public Transform
 BR_REGISTER(Transform, QuantizeTransform)
 
 /*!
+ * \ingroup distances
+ * \brief Bayesian quantization distance
+ * \author Josh Klontz \cite jklontz
+ */
+class BayesianQuantizationDistance : public Distance
+{
+    Q_OBJECT
+    QVector<float> loglikelihood;
+
+    void train(const TemplateList &src)
+    {
+        if (src.first().size() > 1)
+            qFatal("Expected sigle matrix templates.");
+
+        Mat data = OpenCVUtils::toMat(src.data());
+        QList<int> labels = src.labels<int>();
+
+        QVector<qint64> genuines(256*256,0), impostors(256*256,0);
+        for (int i=0; i<labels.size(); i++) {
+            const uchar *a = data.ptr(i);
+            for (int j=0; j<labels.size(); j++) {
+                const uchar *b = data.ptr(j);
+                const bool genuine = (labels[i] == labels[j]);
+                for (int k=0; k<data.cols; k++)
+                    genuine ? genuines[256*a[k]+b[k]]++ : impostors[256*a[k]+b[k]]++;
+            }
+        }
+
+        qint64 totalGenuines(0), totalImpostors(0);
+        for (int i=0; i<256*256; i++) {
+            totalGenuines += genuines[i];
+            totalImpostors += impostors[i];
+        }
+
+        loglikelihood = QVector<float>(256*256);
+        for (int i=0; i<256; i++)
+            for (int j=0; j<256; j++)
+                loglikelihood[i*256+j] = log((double(genuines[i*256+j]+genuines[j*256+i]+1)/totalGenuines)/
+                                             (double(impostors[i*256+j]+impostors[j*256+i]+1)/totalImpostors));
+    }
+
+    float compare(const Template &a, const Template &b) const
+    {
+        const uchar *aData = a.m().data;
+        const uchar *bData = b.m().data;
+        const int size = a.m().rows * a.m().cols;
+        float likelihood = 0;
+        for (int i=0; i<size; i++)
+            likelihood += loglikelihood[256*aData[i]+bData[i]];
+        return likelihood;
+    }
+
+    void load(QDataStream &stream)
+    {
+        stream >> loglikelihood;
+    }
+
+    void store(QDataStream &stream) const
+    {
+        stream << loglikelihood;
+    }
+};
+
+BR_REGISTER(Distance, BayesianQuantizationDistance)
+
+/*!
  * \ingroup transforms
  * \brief Approximate floats as signed bit.
  * \author Josh Klontz \cite jklontz
@@ -169,6 +235,19 @@ public:
     }
 
 private:
+    static void getScores(const QList<int> &indicies, const QList<int> &labels, const Mat &lut, QVector<float> &genuineScores, QVector<float> &impostorScores)
+    {
+        genuineScores.clear(); impostorScores.clear();
+        genuineScores.reserve(indicies.size());
+        impostorScores.reserve(indicies.size()*indicies.size()/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);
+            }
+    }
+
     void _train(const Mat &data, const QList<int> &labels, Mat *lut, Mat *center)
     {
         Mat clusterLabels;
@@ -181,14 +260,17 @@ private:
         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);
-            }
+        QVector<float> genuineScores, impostorScores;
+
+        // RBF Kernel
+//        getScores(indicies, labels, *lut, genuineScores, impostorScores);
+//        float sigma = 1.0 / Common::Mean(impostorScores);
+//        for (int j=0; j<256; j++)
+//            for (int k=0; k<256; k++)
+//                lut->at<float>(0,j*256+k) = exp(-lut->at<float>(0,j*256+k)/(2*pow(sigma, 2.f)));
+
+        // Bayesian PDF
+        getScores(indicies, labels, *lut, genuineScores, impostorScores);
         genuineScores = Common::Downsample(genuineScores, 256);
         impostorScores = Common::Downsample(impostorScores, 256);
 
@@ -199,21 +281,45 @@ private:
             for (int k=0; k<256; k++)
                 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));
+//                lut->at<float>(0,j*256+k) = std::max(0.0, log(Common::KernelDensityEstimation(genuineScores, lut->at<float>(0,j*256+k), hGenuine) /
+//                                                              Common::KernelDensityEstimation(impostorScores, lut->at<float>(0,j*256+k), hImpostor)));
+    }
+
+    int getStep(int cols) const
+    {
+        if (n > 0) return n;
+        if (n == 0) return cols;
+        return ceil(float(cols)/abs(n));
+    }
+
+    int getOffset(int cols) const
+    {
+        if (n >= 0) return 0;
+        const int step = getStep(cols);
+        return (step - cols%step) % step;
+    }
+
+    int getDims(int cols) const
+    {
+        const int step = getStep(cols);
+        if (n >= 0) return cols/step;
+        return ceil(float(cols)/step);
     }
 
     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 int step = getStep(data.cols);
         const QList<int> labels = src.labels<int>();
 
         Mat &lut = ProductQuantizationLUTs[index];
-        lut = Mat(data.cols/n, 256*256, CV_32FC1);
+        lut = Mat(getDims(data.cols), 256*256, CV_32FC1);
 
         QList<Mat> subdata, subluts;
+        const int offset = getOffset(data.cols);
         for (int i=0; i<lut.rows; i++) {
             centers.append(Mat());
-            subdata.append(data.colRange(i*n,(i+1)*n));
+            subdata.append(data.colRange(max(0, i*step-offset), (i+1)*step-offset));
             subluts.append(lut.row(i));
         }
 
@@ -242,9 +348,11 @@ private:
     void project(const Template &src, Template &dst) const
     {
         Mat m = src.m().reshape(1, 1);
-        dst = Mat(1, m.cols/n, CV_8UC1);
+        const int step = getStep(m.cols);
+        const int offset = getOffset(m.cols);
+        dst = Mat(1, getDims(m.cols), CV_8UC1);
         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]);
+            dst.m().at<uchar>(0,i) = getIndex(m.colRange(max(0, i*step-offset), (i+1)*step-offset), centers[i]);
     }
 
     void store(QDataStream &stream) const