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/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/quality.cpp b/openbr/plugins/quality.cpp
index 3187e95..4ff8229 100644
--- a/openbr/plugins/quality.cpp
+++ b/openbr/plugins/quality.cpp
@@ -32,8 +32,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);
     }
 
diff --git a/openbr/plugins/quantize.cpp b/openbr/plugins/quantize.cpp
index 8f0b3cc..5493326 100644
--- a/openbr/plugins/quantize.cpp
+++ b/openbr/plugins/quantize.cpp
@@ -169,6 +169,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 +194,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,6 +215,8 @@ 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)));
     }
 
     void train(const TemplateList &src)