Peter M. Groen / openbr

Commit 7f8093c9934ddc1f73e86bbef3eb62a45ca2ee7a

Authored by DepthDeluxe
1 parent 3505e097

added CUDAPCATransform plugin sketch 

currently uses CPU code from plugins/core but will start to change things one by one
Inline Side-by-side
Showing 1 changed file with 205 additions and 0 deletions
openbr/plugins/cuda/cudapca.cpp 0 → 100644
  View file @7f8093c
  1 +/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
  2 + * Copyright 2012 The MITRE Corporation *
  3 + * *
  4 + * Licensed under the Apache License, Version 2.0 (the "License"); *
  5 + * you may not use this file except in compliance with the License. *
  6 + * You may obtain a copy of the License at *
  7 + * *
  8 + * http://www.apache.org/licenses/LICENSE-2.0 *
  9 + * *
  10 + * Unless required by applicable law or agreed to in writing, software *
  11 + * distributed under the License is distributed on an "AS IS" BASIS, *
  12 + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. *
  13 + * See the License for the specific language governing permissions and *
  14 + * limitations under the License. *
  15 + * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
  16 +// CUDA includes
  17 +#include <cusolverDn.h>
  18 +
  19 +#include <Eigen/Dense>
  20 +#include <openbr/plugins/openbr_internal.h>
  21 +
  22 +#include <openbr/core/common.h>
  23 +#include <openbr/core/eigenutils.h>
  24 +#include <openbr/core/opencvutils.h>
  25 +
  26 +namespace br
  27 +{
  28 +/*!
  29 + * \ingroup transforms
  30 + * \brief Projects input into learned Principal Component Analysis subspace using CUDA.
  31 + * \author Brendan Klare \cite bklare
  32 + * \author Josh Klontz \cite jklontz
  33 + * \author Colin Heinzmann \cite DepthDeluxe
  34 + *
  35 + * \br_property float keep Options are: [keep < 0 - All eigenvalues are retained, keep == 0 - No PCA is performed and the eigenvectors form an identity matrix, 0 < keep < 1 - Keep is the fraction of the variance to retain, keep >= 1 - keep is the number of leading eigenvectors to retain] Default is 0.95.
  36 + * \br_property int drop The number of leading eigen-dimensions to drop.
  37 + * \br_property bool whiten Whether or not to perform PCA whitening (i.e., normalize variance of each dimension to unit norm)
  38 + */
  39 +class CUDAPCATransform : public Transform
  40 +{
  41 + Q_OBJECT
  42 +
  43 +protected:
  44 + Q_PROPERTY(float keep READ get_keep WRITE set_keep RESET reset_keep STORED false)
  45 + Q_PROPERTY(int drop READ get_drop WRITE set_drop RESET reset_drop STORED false)
  46 + Q_PROPERTY(bool whiten READ get_whiten WRITE set_whiten RESET reset_whiten STORED false)
  47 +
  48 + BR_PROPERTY(float, keep, 0.95)
  49 + BR_PROPERTY(int, drop, 0)
  50 + BR_PROPERTY(bool, whiten, false)
  51 +
  52 + Eigen::VectorXf mean, eVals;
  53 + Eigen::MatrixXf eVecs;
  54 +
  55 + int originalRows;
  56 +
  57 +public:
  58 + CUDAPCATransform() : keep(0.95), drop(0), whiten(false) {}
  59 +
  60 +private:
  61 + double residualReconstructionError(const Template &src) const
  62 + {
  63 + Template proj;
  64 + project(src, proj);
  65 +
  66 + Eigen::Map<const Eigen::VectorXf> srcMap(src.m().ptr<float>(), src.m().rows*src.m().cols);
  67 + Eigen::Map<Eigen::VectorXf> projMap(proj.m().ptr<float>(), keep);
  68 +
  69 + return (srcMap - mean).squaredNorm() - projMap.squaredNorm();
  70 + }
  71 +
  72 + void train(const TemplateList &trainingSet)
  73 + {
  74 + if (trainingSet.first().m().type() != CV_32FC1)
  75 + qFatal("Requires single channel 32-bit floating point matrices.");
  76 +
  77 + originalRows = trainingSet.first().m().rows;
  78 + int dimsIn = trainingSet.first().m().rows * trainingSet.first().m().cols;
  79 + const int instances = trainingSet.size();
  80 +
  81 + // Map into 64-bit Eigen matrix
  82 + Eigen::MatrixXd data(dimsIn, instances);
  83 + for (int i=0; i<instances; i++)
  84 + data.col(i) = Eigen::Map<const Eigen::MatrixXf>(trainingSet[i].m().ptr<float>(), dimsIn, 1).cast<double>();
  85 +
  86 + trainCore(data);
  87 + }
  88 +
  89 + void project(const Template &src, Template &dst) const
  90 + {
  91 + dst = cv::Mat(1, keep, CV_32FC1);
  92 +
  93 + // Map Eigen into OpenCV
  94 + Eigen::Map<const Eigen::MatrixXf> inMap(src.m().ptr<float>(), src.m().rows*src.m().cols, 1);
  95 + Eigen::Map<Eigen::MatrixXf> outMap(dst.m().ptr<float>(), keep, 1);
  96 +
  97 + // Do projection
  98 + outMap = eVecs.transpose() * (inMap - mean);
  99 + }
  100 +
  101 + void store(QDataStream &stream) const
  102 + {
  103 + stream << keep << drop << whiten << originalRows << mean << eVals << eVecs;
  104 + }
  105 +
  106 + void load(QDataStream &stream)
  107 + {
  108 + stream >> keep >> drop >> whiten >> originalRows >> mean >> eVals >> eVecs;
  109 + }
  110 +
  111 +protected:
  112 + void trainCore(Eigen::MatrixXd data)
  113 + {
  114 + int dimsIn = data.rows();
  115 + int instances = data.cols();
  116 + const bool dominantEigenEstimation = (dimsIn > instances);
  117 +
  118 + Eigen::MatrixXd allEVals, allEVecs;
  119 + if (keep != 0) {
  120 + // Compute and remove mean
  121 + mean = Eigen::VectorXf(dimsIn);
  122 + for (int i=0; i<dimsIn; i++) mean(i) = data.row(i).sum() / (float)instances;
  123 + for (int i=0; i<dimsIn; i++) data.row(i).array() -= mean(i);
  124 +
  125 + // Calculate covariance matrix
  126 + Eigen::MatrixXd cov;
  127 + if (dominantEigenEstimation) cov = data.transpose() * data / (instances-1.0);
  128 + else cov = data * data.transpose() / (instances-1.0);
  129 +
  130 + // Compute eigendecomposition. Returns eigenvectors/eigenvalues in increasing order by eigenvalue.
  131 + Eigen::SelfAdjointEigenSolver<Eigen::MatrixXd> eSolver(cov);
  132 + allEVals = eSolver.eigenvalues();
  133 + allEVecs = eSolver.eigenvectors();
  134 + if (dominantEigenEstimation) allEVecs = data * allEVecs;
  135 + } else {
  136 + // Null case
  137 + mean = Eigen::VectorXf::Zero(dimsIn);
  138 + allEVecs = Eigen::MatrixXd::Identity(dimsIn, dimsIn);
  139 + allEVals = Eigen::VectorXd::Ones(dimsIn);
  140 + }
  141 +
  142 + if (keep <= 0) {
  143 + keep = dimsIn - drop;
  144 + } else if (keep < 1) {
  145 + // Keep eigenvectors that retain a certain energy percentage.
  146 + const double totalEnergy = allEVals.sum();
  147 + if (totalEnergy == 0) {
  148 + keep = 0;
  149 + } else {
  150 + double currentEnergy = 0;
  151 + int i=0;
  152 + while ((currentEnergy / totalEnergy < keep) && (i < allEVals.rows())) {
  153 + currentEnergy += allEVals(allEVals.rows()-(i+1));
  154 + i++;
  155 + }
  156 + keep = i - drop;
  157 + }
  158 + } else {
  159 + if (keep + drop > allEVals.rows()) {
  160 + qWarning("Insufficient samples, needed at least %d but only got %d.", (int)keep + drop, (int)allEVals.rows());
  161 + keep = allEVals.rows() - drop;
  162 + }
  163 + }
  164 +
  165 + // Keep highest energy vectors
  166 + eVals = Eigen::VectorXf((int)keep, 1);
  167 + eVecs = Eigen::MatrixXf(allEVecs.rows(), (int)keep);
  168 + for (int i=0; i<keep; i++) {
  169 + int index = allEVals.rows()-(i+drop+1);
  170 + eVals(i) = allEVals(index);
  171 + eVecs.col(i) = allEVecs.col(index).cast<float>() / allEVecs.col(index).norm();
  172 + if (whiten) eVecs.col(i) /= sqrt(eVals(i));
  173 + }
  174 +
  175 + // Debug output
  176 + if (Globals->verbose) qDebug() << "PCA Training:\n\tDimsIn =" << dimsIn << "\n\tKeep =" << keep;
  177 + }
  178 +
  179 + void writeEigenVectors(const Eigen::MatrixXd &allEVals, const Eigen::MatrixXd &allEVecs) const
  180 + {
  181 + const int originalCols = mean.rows() / originalRows;
  182 +
  183 + { // Write out mean image
  184 + cv::Mat out(originalRows, originalCols, CV_32FC1);
  185 + Eigen::Map<Eigen::MatrixXf> outMap(out.ptr<float>(), mean.rows(), 1);
  186 + outMap = mean.col(0);
  187 + // OpenCVUtils::saveImage(out, Globals->Debug+"/PCA/eigenVectors/mean.png");
  188 + }
  189 +
  190 + // Write out sample eigen vectors (16 highest, 8 lowest), filename = eigenvalue.
  191 + for (int k=0; k<(int)allEVals.size(); k++) {
  192 + if ((k < 8) || (k >= (int)allEVals.size()-16)) {
  193 + cv::Mat out(originalRows, originalCols, CV_64FC1);
  194 + Eigen::Map<Eigen::MatrixXd> outMap(out.ptr<double>(), mean.rows(), 1);
  195 + outMap = allEVecs.col(k);
  196 + // OpenCVUtils::saveImage(out, Globals->Debug+"/PCA/eigenVectors/"+QString::number(allEVals(k),'f',0)+".png");
  197 + }
  198 + }
  199 + }
  200 +};
  201 +
  202 +BR_REGISTER(Transform, CUDAPCATransform)
  203 +} // namespace br
  204 +
  205 +#include "cuda/cudapca.moc"