modules/contrib/src/facerec.cpp

   1 /*
   2  * Copyright (c) 2011,2012. Philipp Wagner <bytefish[at]gmx[dot]de>.
   3  * Released to public domain under terms of the BSD Simplified license.
   4  *
   5  * Redistribution and use in source and binary forms, with or without
   6  * modification, are permitted provided that the following conditions are met:
   7  *   * Redistributions of source code must retain the above copyright
   8  *     notice, this list of conditions and the following disclaimer.
   9  *   * Redistributions in binary form must reproduce the above copyright
  10  *     notice, this list of conditions and the following disclaimer in the
  11  *     documentation and/or other materials provided with the distribution.
  12  *   * Neither the name of the organization nor the names of its contributors
  13  *     may be used to endorse or promote products derived from this software
  14  *     without specific prior written permission.
  15  *
  16  *   See <http://www.opensource.org/licenses/bsd-license>
  17  */
  18 #include "precomp.hpp"
  19 #include <set>
  20 #include <limits>
  21
  22 namespace cv
  23 {
  24
  25 // Reads a sequence from a FileNode::SEQ with type _Tp into a result vector.
  26 template<typename _Tp>
  27 inline void readFileNodeList(const FileNode& fn, std::vector<_Tp>& result) {
  28     if (fn.type() == FileNode::SEQ) {
  29         for (FileNodeIterator it = fn.begin(); it != fn.end();) {
  30             _Tp item;
  31             it >> item;
  32             result.push_back(item);
  33         }
  34     }
  35 }
  36
  37 // Writes the a list of given items to a cv::FileStorage.
  38 template<typename _Tp>
  39 inline void writeFileNodeList(FileStorage& fs, const String& name,
  40                               const std::vector<_Tp>& items) {
  41     // typedefs
  42     typedef typename std::vector<_Tp>::const_iterator constVecIterator;
  43     // write the elements in item to fs
  44     fs << name << "[";
  45     for (constVecIterator it = items.begin(); it != items.end(); ++it) {
  46         fs << *it;
  47     }
  48     fs << "]";
  49 }
  50
  51 static Mat asRowMatrix(InputArrayOfArrays src, int rtype, double alpha=1, double beta=0) {
  52     // make sure the input data is a vector of matrices or vector of vector
  53     if(src.kind() != _InputArray::STD_VECTOR_MAT && src.kind() != _InputArray::STD_VECTOR_VECTOR) {
  54         String error_message = "The data is expected as InputArray::STD_VECTOR_MAT (a std::vector<Mat>) or _InputArray::STD_VECTOR_VECTOR (a std::vector< std::vector<...> >).";
  55         CV_Error(Error::StsBadArg, error_message);
  56     }
  57     // number of samples
  58     size_t n = src.total();
  59     // return empty matrix if no matrices given
  60     if(n == 0)
  61         return Mat();
  62     // dimensionality of (reshaped) samples
  63     size_t d = src.getMat(0).total();
  64     // create data matrix
  65     Mat data((int)n, (int)d, rtype);
  66     // now copy data
  67     for(unsigned int i = 0; i < n; i++) {
  68         // make sure data can be reshaped, throw exception if not!
  69         if(src.getMat(i).total() != d) {
  70             String error_message = format("Wrong number of elements in matrix #%d! Expected %d was %d.", i, d, src.getMat(i).total());
  71             CV_Error(Error::StsBadArg, error_message);
  72         }
  73         // get a hold of the current row
  74         Mat xi = data.row(i);
  75         // make reshape happy by cloning for non-continuous matrices
  76         if(src.getMat(i).isContinuous()) {
  77             src.getMat(i).reshape(1, 1).convertTo(xi, rtype, alpha, beta);
  78         } else {
  79             src.getMat(i).clone().reshape(1, 1).convertTo(xi, rtype, alpha, beta);
  80         }
  81     }
  82     return data;
  83 }
  84
  85
  86 // Removes duplicate elements in a given vector.
  87 template<typename _Tp>
  88 inline std::vector<_Tp> remove_dups(const std::vector<_Tp>& src) {
  89     typedef typename std::set<_Tp>::const_iterator constSetIterator;
  90     typedef typename std::vector<_Tp>::const_iterator constVecIterator;
  91     std::set<_Tp> set_elems;
  92     for (constVecIterator it = src.begin(); it != src.end(); ++it)
  93         set_elems.insert(*it);
  94     std::vector<_Tp> elems;
  95     for (constSetIterator it = set_elems.begin(); it != set_elems.end(); ++it)
  96         elems.push_back(*it);
  97     return elems;
  98 }
  99
 100
 101 // Turk, M., and Pentland, A. "Eigenfaces for recognition.". Journal of
 102 // Cognitive Neuroscience 3 (1991), 71–86.
 103 class Eigenfaces : public FaceRecognizer
 104 {
 105 private:
 106     int _num_components;
 107     double _threshold;
 108     std::vector<Mat> _projections;
 109     Mat _labels;
 110     Mat _eigenvectors;
 111     Mat _eigenvalues;
 112     Mat _mean;
 113
 114 public:
 115     using FaceRecognizer::save;
 116     using FaceRecognizer::load;
 117
 118     // Initializes an empty Eigenfaces model.
 119     Eigenfaces(int num_components = 0, double threshold = DBL_MAX) :
 120         _num_components(num_components),
 121         _threshold(threshold) {}
 122
 123     // Initializes and computes an Eigenfaces model with images in src and
 124     // corresponding labels in labels. num_components will be kept for
 125     // classification.
 126     Eigenfaces(InputArrayOfArrays src, InputArray labels,
 127             int num_components = 0, double threshold = DBL_MAX) :
 128         _num_components(num_components),
 129         _threshold(threshold) {
 130         train(src, labels);
 131     }
 132
 133     // Computes an Eigenfaces model with images in src and corresponding labels
 134     // in labels.
 135     void train(InputArrayOfArrays src, InputArray labels);
 136
 137     // Predicts the label of a query image in src.
 138     int predict(InputArray src) const;
 139
 140     // Predicts the label and confidence for a given sample.
 141     void predict(InputArray _src, int &label, double &dist) const;
 142
 143     // See FaceRecognizer::load.
 144     void load(const FileStorage& fs);
 145
 146     // See FaceRecognizer::save.
 147     void save(FileStorage& fs) const;
 148
 149     AlgorithmInfo* info() const;
 150 };
 151
 152 // Belhumeur, P. N., Hespanha, J., and Kriegman, D. "Eigenfaces vs. Fisher-
 153 // faces: Recognition using class specific linear projection.". IEEE
 154 // Transactions on Pattern Analysis and Machine Intelligence 19, 7 (1997),
 155 // 711–720.
 156 class Fisherfaces: public FaceRecognizer
 157 {
 158 private:
 159     int _num_components;
 160     double _threshold;
 161     Mat _eigenvectors;
 162     Mat _eigenvalues;
 163     Mat _mean;
 164     std::vector<Mat> _projections;
 165     Mat _labels;
 166
 167 public:
 168     using FaceRecognizer::save;
 169     using FaceRecognizer::load;
 170
 171     // Initializes an empty Fisherfaces model.
 172     Fisherfaces(int num_components = 0, double threshold = DBL_MAX) :
 173         _num_components(num_components),
 174         _threshold(threshold) {}
 175
 176     // Initializes and computes a Fisherfaces model with images in src and
 177     // corresponding labels in labels. num_components will be kept for
 178     // classification.
 179     Fisherfaces(InputArrayOfArrays src, InputArray labels,
 180             int num_components = 0, double threshold = DBL_MAX) :
 181         _num_components(num_components),
 182         _threshold(threshold) {
 183         train(src, labels);
 184     }
 185
 186     ~Fisherfaces() {}
 187
 188     // Computes a Fisherfaces model with images in src and corresponding labels
 189     // in labels.
 190     void train(InputArrayOfArrays src, InputArray labels);
 191
 192     // Predicts the label of a query image in src.
 193     int predict(InputArray src) const;
 194
 195     // Predicts the label and confidence for a given sample.
 196     void predict(InputArray _src, int &label, double &dist) const;
 197
 198     // See FaceRecognizer::load.
 199     void load(const FileStorage& fs);
 200
 201     // See FaceRecognizer::save.
 202     void save(FileStorage& fs) const;
 203
 204     AlgorithmInfo* info() const;
 205 };
 206
 207 // Face Recognition based on Local Binary Patterns.
 208 //
 209 //  Ahonen T, Hadid A. and Pietikäinen M. "Face description with local binary
 210 //  patterns: Application to face recognition." IEEE Transactions on Pattern
 211 //  Analysis and Machine Intelligence, 28(12):2037-2041.
 212 //
 213 class LBPH : public FaceRecognizer
 214 {
 215 private:
 216     int _grid_x;
 217     int _grid_y;
 218     int _radius;
 219     int _neighbors;
 220     double _threshold;
 221
 222     std::vector<Mat> _histograms;
 223     Mat _labels;
 224
 225     // Computes a LBPH model with images in src and
 226     // corresponding labels in labels, possibly preserving
 227     // old model data.
 228     void train(InputArrayOfArrays src, InputArray labels, bool preserveData);
 229
 230
 231 public:
 232     using FaceRecognizer::save;
 233     using FaceRecognizer::load;
 234
 235     // Initializes this LBPH Model. The current implementation is rather fixed
 236     // as it uses the Extended Local Binary Patterns per default.
 237     //
 238     // radius, neighbors are used in the local binary patterns creation.
 239     // grid_x, grid_y control the grid size of the spatial histograms.
 240     LBPH(int radius_=1, int neighbors_=8,
 241             int gridx=8, int gridy=8,
 242             double threshold = DBL_MAX) :
 243         _grid_x(gridx),
 244         _grid_y(gridy),
 245         _radius(radius_),
 246         _neighbors(neighbors_),
 247         _threshold(threshold) {}
 248
 249     // Initializes and computes this LBPH Model. The current implementation is
 250     // rather fixed as it uses the Extended Local Binary Patterns per default.
 251     //
 252     // (radius=1), (neighbors=8) are used in the local binary patterns creation.
 253     // (grid_x=8), (grid_y=8) controls the grid size of the spatial histograms.
 254     LBPH(InputArrayOfArrays src,
 255             InputArray labels,
 256             int radius_=1, int neighbors_=8,
 257             int gridx=8, int gridy=8,
 258             double threshold = DBL_MAX) :
 259                 _grid_x(gridx),
 260                 _grid_y(gridy),
 261                 _radius(radius_),
 262                 _neighbors(neighbors_),
 263                 _threshold(threshold) {
 264         train(src, labels);
 265     }
 266
 267     ~LBPH() { }
 268
 269     // Computes a LBPH model with images in src and
 270     // corresponding labels in labels.
 271     void train(InputArrayOfArrays src, InputArray labels);
 272
 273     // Updates this LBPH model with images in src and
 274     // corresponding labels in labels.
 275     void update(InputArrayOfArrays src, InputArray labels);
 276
 277     // Predicts the label of a query image in src.
 278     int predict(InputArray src) const;
 279
 280     // Predicts the label and confidence for a given sample.
 281     void predict(InputArray _src, int &label, double &dist) const;
 282
 283     // See FaceRecognizer::load.
 284     void load(const FileStorage& fs);
 285
 286     // See FaceRecognizer::save.
 287     void save(FileStorage& fs) const;
 288
 289     // Getter functions.
 290     int neighbors() const { return _neighbors; }
 291     int radius() const { return _radius; }
 292     int grid_x() const { return _grid_x; }
 293     int grid_y() const { return _grid_y; }
 294
 295     AlgorithmInfo* info() const;
 296 };
 297
 298
 299 //------------------------------------------------------------------------------
 300 // FaceRecognizer
 301 //------------------------------------------------------------------------------
 302 void FaceRecognizer::update(InputArrayOfArrays src, InputArray labels ) {
 303     if( dynamic_cast<LBPH*>(this) != 0 )
 304     {
 305         dynamic_cast<LBPH*>(this)->update( src, labels );
 306         return;
 307     }
 308
 309     String error_msg = format("This FaceRecognizer (%s) does not support updating, you have to use FaceRecognizer::train to update it.", this->name().c_str());
 310     CV_Error(Error::StsNotImplemented, error_msg);
 311 }
 312
 313 void FaceRecognizer::save(const String& filename) const {
 314     FileStorage fs(filename, FileStorage::WRITE);
 315     if (!fs.isOpened())
 316         CV_Error(Error::StsError, "File can't be opened for writing!");
 317     this->save(fs);
 318     fs.release();
 319 }
 320
 321 void FaceRecognizer::load(const String& filename) {
 322     FileStorage fs(filename, FileStorage::READ);
 323     if (!fs.isOpened())
 324         CV_Error(Error::StsError, "File can't be opened for writing!");
 325     this->load(fs);
 326     fs.release();
 327 }
 328
 329 //------------------------------------------------------------------------------
 330 // Eigenfaces
 331 //------------------------------------------------------------------------------
 332 void Eigenfaces::train(InputArrayOfArrays _src, InputArray _local_labels) {
 333     if(_src.total() == 0) {
 334         String error_message = format("Empty training data was given. You'll need more than one sample to learn a model.");
 335         CV_Error(Error::StsBadArg, error_message);
 336     } else if(_local_labels.getMat().type() != CV_32SC1) {
 337         String error_message = format("Labels must be given as integer (CV_32SC1). Expected %d, but was %d.", CV_32SC1, _local_labels.type());
 338         CV_Error(Error::StsBadArg, error_message);
 339     }
 340     // make sure data has correct size
 341     if(_src.total() > 1) {
 342         for(int i = 1; i < static_cast<int>(_src.total()); i++) {
 343             if(_src.getMat(i-1).total() != _src.getMat(i).total()) {
 344                 String error_message = format("In the Eigenfaces method all input samples (training images) must be of equal size! Expected %d pixels, but was %d pixels.", _src.getMat(i-1).total(), _src.getMat(i).total());
 345                 CV_Error(Error::StsUnsupportedFormat, error_message);
 346             }
 347         }
 348     }
 349     // get labels
 350     Mat labels = _local_labels.getMat();
 351     // observations in row
 352     Mat data = asRowMatrix(_src, CV_64FC1);
 353
 354     // number of samples
 355    int n = data.rows;
 356     // assert there are as much samples as labels
 357     if(static_cast<int>(labels.total()) != n) {
 358         String error_message = format("The number of samples (src) must equal the number of labels (labels)! len(src)=%d, len(labels)=%d.", n, labels.total());
 359         CV_Error(Error::StsBadArg, error_message);
 360     }
 361     // clear existing model data
 362     _labels.release();
 363     _projections.clear();
 364     // clip number of components to be valid
 365     if((_num_components <= 0) || (_num_components > n))
 366         _num_components = n;
 367
 368     // perform the PCA
 369     PCA pca(data, Mat(), PCA::DATA_AS_ROW, _num_components);
 370     // copy the PCA results
 371     _mean = pca.mean.reshape(1,1); // store the mean vector
 372     _eigenvalues = pca.eigenvalues.clone(); // eigenvalues by row
 373     transpose(pca.eigenvectors, _eigenvectors); // eigenvectors by column
 374     // store labels for prediction
 375     _labels = labels.clone();
 376     // save projections
 377     for(int sampleIdx = 0; sampleIdx < data.rows; sampleIdx++) {
 378         Mat p = subspaceProject(_eigenvectors, _mean, data.row(sampleIdx));
 379         _projections.push_back(p);
 380     }
 381 }
 382
 383 void Eigenfaces::predict(InputArray _src, int &minClass, double &minDist) const {
 384     // get data
 385     Mat src = _src.getMat();
 386     // make sure the user is passing correct data
 387     if(_projections.empty()) {
 388         // throw error if no data (or simply return -1?)
 389         String error_message = "This Eigenfaces model is not computed yet. Did you call Eigenfaces::train?";
 390         CV_Error(Error::StsError, error_message);
 391     } else if(_eigenvectors.rows != static_cast<int>(src.total())) {
 392         // check data alignment just for clearer exception messages
 393         String error_message = format("Wrong input image size. Reason: Training and Test images must be of equal size! Expected an image with %d elements, but got %d.", _eigenvectors.rows, src.total());
 394         CV_Error(Error::StsBadArg, error_message);
 395     }
 396     // project into PCA subspace
 397     Mat q = subspaceProject(_eigenvectors, _mean, src.reshape(1,1));
 398     minDist = DBL_MAX;
 399     minClass = -1;
 400     for(size_t sampleIdx = 0; sampleIdx < _projections.size(); sampleIdx++) {
 401         double dist = norm(_projections[sampleIdx], q, NORM_L2);
 402         if((dist < minDist) && (dist < _threshold)) {
 403             minDist = dist;
 404             minClass = _labels.at<int>((int)sampleIdx);
 405         }
 406     }
 407 }
 408
 409 int Eigenfaces::predict(InputArray _src) const {
 410     int label;
 411     double dummy;
 412     predict(_src, label, dummy);
 413     return label;
 414 }
 415
 416 void Eigenfaces::load(const FileStorage& fs) {
 417     //read matrices
 418     fs["num_components"] >> _num_components;
 419     fs["mean"] >> _mean;
 420     fs["eigenvalues"] >> _eigenvalues;
 421     fs["eigenvectors"] >> _eigenvectors;
 422     // read sequences
 423     readFileNodeList(fs["projections"], _projections);
 424     fs["labels"] >> _labels;
 425 }
 426
 427 void Eigenfaces::save(FileStorage& fs) const {
 428     // write matrices
 429     fs << "num_components" << _num_components;
 430     fs << "mean" << _mean;
 431     fs << "eigenvalues" << _eigenvalues;
 432     fs << "eigenvectors" << _eigenvectors;
 433     // write sequences
 434     writeFileNodeList(fs, "projections", _projections);
 435     fs << "labels" << _labels;
 436 }
 437
 438 //------------------------------------------------------------------------------
 439 // Fisherfaces
 440 //------------------------------------------------------------------------------
 441 void Fisherfaces::train(InputArrayOfArrays src, InputArray _lbls) {
 442     if(src.total() == 0) {
 443         String error_message = format("Empty training data was given. You'll need more than one sample to learn a model.");
 444         CV_Error(Error::StsBadArg, error_message);
 445     } else if(_lbls.getMat().type() != CV_32SC1) {
 446         String error_message = format("Labels must be given as integer (CV_32SC1). Expected %d, but was %d.", CV_32SC1, _lbls.type());
 447         CV_Error(Error::StsBadArg, error_message);
 448     }
 449     // make sure data has correct size
 450     if(src.total() > 1) {
 451         for(int i = 1; i < static_cast<int>(src.total()); i++) {
 452             if(src.getMat(i-1).total() != src.getMat(i).total()) {
 453                 String error_message = format("In the Fisherfaces method all input samples (training images) must be of equal size! Expected %d pixels, but was %d pixels.", src.getMat(i-1).total(), src.getMat(i).total());
 454                 CV_Error(Error::StsUnsupportedFormat, error_message);
 455             }
 456         }
 457     }
 458     // get data
 459     Mat labels = _lbls.getMat();
 460     Mat data = asRowMatrix(src, CV_64FC1);
 461     // number of samples
 462     int N = data.rows;
 463     // make sure labels are passed in correct shape
 464     if(labels.total() != (size_t) N) {
 465         String error_message = format("The number of samples (src) must equal the number of labels (labels)! len(src)=%d, len(labels)=%d.", N, labels.total());
 466         CV_Error(Error::StsBadArg, error_message);
 467     } else if(labels.rows != 1 && labels.cols != 1) {
 468         String error_message = format("Expected the labels in a matrix with one row or column! Given dimensions are rows=%s, cols=%d.", labels.rows, labels.cols);
 469        CV_Error(Error::StsBadArg, error_message);
 470     }
 471     // clear existing model data
 472     _labels.release();
 473     _projections.clear();
 474     // safely copy from cv::Mat to std::vector
 475     std::vector<int> ll;
 476     for(unsigned int i = 0; i < labels.total(); i++) {
 477         ll.push_back(labels.at<int>(i));
 478     }
 479     // get the number of unique classes
 480     int C = (int) remove_dups(ll).size();
 481     // clip number of components to be a valid number
 482     if((_num_components <= 0) || (_num_components > (C-1)))
 483         _num_components = (C-1);
 484     // perform a PCA and keep (N-C) components
 485     PCA pca(data, Mat(), PCA::DATA_AS_ROW, (N-C));
 486     // project the data and perform a LDA on it
 487     LDA lda(pca.project(data),labels, _num_components);
 488     // store the total mean vector
 489     _mean = pca.mean.reshape(1,1);
 490     // store labels
 491     _labels = labels.clone();
 492     // store the eigenvalues of the discriminants
 493     lda.eigenvalues().convertTo(_eigenvalues, CV_64FC1);
 494     // Now calculate the projection matrix as pca.eigenvectors * lda.eigenvectors.
 495     // Note: OpenCV stores the eigenvectors by row, so we need to transpose it!
 496     gemm(pca.eigenvectors, lda.eigenvectors(), 1.0, Mat(), 0.0, _eigenvectors, GEMM_1_T);
 497     // store the projections of the original data
 498     for(int sampleIdx = 0; sampleIdx < data.rows; sampleIdx++) {
 499         Mat p = subspaceProject(_eigenvectors, _mean, data.row(sampleIdx));
 500         _projections.push_back(p);
 501     }
 502 }
 503
 504 void Fisherfaces::predict(InputArray _src, int &minClass, double &minDist) const {
 505     Mat src = _src.getMat();
 506     // check data alignment just for clearer exception messages
 507     if(_projections.empty()) {
 508         // throw error if no data (or simply return -1?)
 509         String error_message = "This Fisherfaces model is not computed yet. Did you call Fisherfaces::train?";
 510         CV_Error(Error::StsBadArg, error_message);
 511     } else if(src.total() != (size_t) _eigenvectors.rows) {
 512         String error_message = format("Wrong input image size. Reason: Training and Test images must be of equal size! Expected an image with %d elements, but got %d.", _eigenvectors.rows, src.total());
 513         CV_Error(Error::StsBadArg, error_message);
 514     }
 515     // project into LDA subspace
 516     Mat q = subspaceProject(_eigenvectors, _mean, src.reshape(1,1));
 517     // find 1-nearest neighbor
 518     minDist = DBL_MAX;
 519     minClass = -1;
 520     for(size_t sampleIdx = 0; sampleIdx < _projections.size(); sampleIdx++) {
 521         double dist = norm(_projections[sampleIdx], q, NORM_L2);
 522         if((dist < minDist) && (dist < _threshold)) {
 523             minDist = dist;
 524             minClass = _labels.at<int>((int)sampleIdx);
 525         }
 526     }
 527 }
 528
 529 int Fisherfaces::predict(InputArray _src) const {
 530     int label;
 531     double dummy;
 532     predict(_src, label, dummy);
 533     return label;
 534 }
 535
 536 // See FaceRecognizer::load.
 537 void Fisherfaces::load(const FileStorage& fs) {
 538     //read matrices
 539     fs["num_components"] >> _num_components;
 540     fs["mean"] >> _mean;
 541     fs["eigenvalues"] >> _eigenvalues;
 542     fs["eigenvectors"] >> _eigenvectors;
 543     // read sequences
 544     readFileNodeList(fs["projections"], _projections);
 545     fs["labels"] >> _labels;
 546 }
 547
 548 // See FaceRecognizer::save.
 549 void Fisherfaces::save(FileStorage& fs) const {
 550     // write matrices
 551     fs << "num_components" << _num_components;
 552     fs << "mean" << _mean;
 553     fs << "eigenvalues" << _eigenvalues;
 554     fs << "eigenvectors" << _eigenvectors;
 555     // write sequences
 556     writeFileNodeList(fs, "projections", _projections);
 557     fs << "labels" << _labels;
 558 }
 559
 560 //------------------------------------------------------------------------------
 561 // LBPH
 562 //------------------------------------------------------------------------------
 563
 564 template <typename _Tp> static
 565 void olbp_(InputArray _src, OutputArray _dst) {
 566     // get matrices
 567     Mat src = _src.getMat();
 568     // allocate memory for result
 569     _dst.create(src.rows-2, src.cols-2, CV_8UC1);
 570     Mat dst = _dst.getMat();
 571     // zero the result matrix
 572     dst.setTo(0);
 573     // calculate patterns
 574     for(int i=1;i<src.rows-1;i++) {
 575         for(int j=1;j<src.cols-1;j++) {
 576             _Tp center = src.at<_Tp>(i,j);
 577             unsigned char code = 0;
 578             code |= (src.at<_Tp>(i-1,j-1) >= center) << 7;
 579             code |= (src.at<_Tp>(i-1,j) >= center) << 6;
 580             code |= (src.at<_Tp>(i-1,j+1) >= center) << 5;
 581             code |= (src.at<_Tp>(i,j+1) >= center) << 4;
 582             code |= (src.at<_Tp>(i+1,j+1) >= center) << 3;
 583             code |= (src.at<_Tp>(i+1,j) >= center) << 2;
 584             code |= (src.at<_Tp>(i+1,j-1) >= center) << 1;
 585             code |= (src.at<_Tp>(i,j-1) >= center) << 0;
 586             dst.at<unsigned char>(i-1,j-1) = code;
 587         }
 588     }
 589 }
 590
 591 //------------------------------------------------------------------------------
 592 // cv::elbp
 593 //------------------------------------------------------------------------------
 594 template <typename _Tp> static
 595 inline void elbp_(InputArray _src, OutputArray _dst, int radius, int neighbors) {
 596     //get matrices
 597     Mat src = _src.getMat();
 598     // allocate memory for result
 599     _dst.create(src.rows-2*radius, src.cols-2*radius, CV_32SC1);
 600     Mat dst = _dst.getMat();
 601     // zero
 602     dst.setTo(0);
 603     for(int n=0; n<neighbors; n++) {
 604         // sample points
 605         float x = static_cast<float>(radius * cos(2.0*CV_PI*n/static_cast<float>(neighbors)));
 606         float y = static_cast<float>(-radius * sin(2.0*CV_PI*n/static_cast<float>(neighbors)));
 607         // relative indices
 608         int fx = static_cast<int>(floor(x));
 609         int fy = static_cast<int>(floor(y));
 610         int cx = static_cast<int>(ceil(x));
 611         int cy = static_cast<int>(ceil(y));
 612         // fractional part
 613         float ty = y - fy;
 614         float tx = x - fx;
 615         // set interpolation weights
 616         float w1 = (1 - tx) * (1 - ty);
 617         float w2 =      tx  * (1 - ty);
 618         float w3 = (1 - tx) *      ty;
 619         float w4 =      tx  *      ty;
 620         // iterate through your data
 621         for(int i=radius; i < src.rows-radius;i++) {
 622             for(int j=radius;j < src.cols-radius;j++) {
 623                 // calculate interpolated value
 624                 float t = static_cast<float>(w1*src.at<_Tp>(i+fy,j+fx) + w2*src.at<_Tp>(i+fy,j+cx) + w3*src.at<_Tp>(i+cy,j+fx) + w4*src.at<_Tp>(i+cy,j+cx));
 625                 // floating point precision, so check some machine-dependent epsilon
 626                 dst.at<int>(i-radius,j-radius) += ((t > src.at<_Tp>(i,j)) || (std::abs(t-src.at<_Tp>(i,j)) < std::numeric_limits<float>::epsilon())) << n;
 627             }
 628         }
 629     }
 630 }
 631
 632 static void elbp(InputArray src, OutputArray dst, int radius, int neighbors)
 633 {
 634     int type = src.type();
 635     switch (type) {
 636     case CV_8SC1:   elbp_<char>(src,dst, radius, neighbors); break;
 637     case CV_8UC1:   elbp_<unsigned char>(src, dst, radius, neighbors); break;
 638     case CV_16SC1:  elbp_<short>(src,dst, radius, neighbors); break;
 639     case CV_16UC1:  elbp_<unsigned short>(src,dst, radius, neighbors); break;
 640     case CV_32SC1:  elbp_<int>(src,dst, radius, neighbors); break;
 641     case CV_32FC1:  elbp_<float>(src,dst, radius, neighbors); break;
 642     case CV_64FC1:  elbp_<double>(src,dst, radius, neighbors); break;
 643     default:
 644         String error_msg = format("Using Original Local Binary Patterns for feature extraction only works on single-channel images (given %d). Please pass the image data as a grayscale image!", type);
 645         CV_Error(Error::StsNotImplemented, error_msg);
 646         break;
 647     }
 648 }
 649
 650 static Mat
 651 histc_(const Mat& src, int minVal=0, int maxVal=255, bool normed=false)
 652 {
 653     Mat result;
 654     // Establish the number of bins.
 655     int histSize = maxVal-minVal+1;
 656     // Set the ranges.
 657     float range[] = { static_cast<float>(minVal), static_cast<float>(maxVal+1) };
 658     const float* histRange = { range };
 659     // calc histogram
 660     calcHist(&src, 1, 0, Mat(), result, 1, &histSize, &histRange, true, false);
 661     // normalize
 662     if(normed) {
 663         result /= (int)src.total();
 664     }
 665     return result.reshape(1,1);
 666 }
 667
 668 static Mat histc(InputArray _src, int minVal, int maxVal, bool normed)
 669 {
 670     Mat src = _src.getMat();
 671     switch (src.type()) {
 672         case CV_8SC1:
 673             return histc_(Mat_<float>(src), minVal, maxVal, normed);
 674             break;
 675         case CV_8UC1:
 676             return histc_(src, minVal, maxVal, normed);
 677             break;
 678         case CV_16SC1:
 679             return histc_(Mat_<float>(src), minVal, maxVal, normed);
 680             break;
 681         case CV_16UC1:
 682             return histc_(src, minVal, maxVal, normed);
 683             break;
 684         case CV_32SC1:
 685             return histc_(Mat_<float>(src), minVal, maxVal, normed);
 686             break;
 687         case CV_32FC1:
 688             return histc_(src, minVal, maxVal, normed);
 689             break;
 690         default:
 691             CV_Error(Error::StsUnmatchedFormats, "This type is not implemented yet."); break;
 692     }
 693     return Mat();
 694 }
 695
 696
 697 static Mat spatial_histogram(InputArray _src, int numPatterns,
 698                              int grid_x, int grid_y, bool /*normed*/)
 699 {
 700     Mat src = _src.getMat();
 701     // calculate LBP patch size
 702     int width = src.cols/grid_x;
 703     int height = src.rows/grid_y;
 704     // allocate memory for the spatial histogram
 705     Mat result = Mat::zeros(grid_x * grid_y, numPatterns, CV_32FC1);
 706     // return matrix with zeros if no data was given
 707     if(src.empty())
 708         return result.reshape(1,1);
 709     // initial result_row
 710     int resultRowIdx = 0;
 711     // iterate through grid
 712     for(int i = 0; i < grid_y; i++) {
 713         for(int j = 0; j < grid_x; j++) {
 714             Mat src_cell = Mat(src, Range(i*height,(i+1)*height), Range(j*width,(j+1)*width));
 715             Mat cell_hist = histc(src_cell, 0, (numPatterns-1), true);
 716             // copy to the result matrix
 717             Mat result_row = result.row(resultRowIdx);
 718             cell_hist.reshape(1,1).convertTo(result_row, CV_32FC1);
 719             // increase row count in result matrix
 720             resultRowIdx++;
 721         }
 722     }
 723     // return result as reshaped feature vector
 724     return result.reshape(1,1);
 725 }
 726
 727 //------------------------------------------------------------------------------
 728 // wrapper to cv::elbp (extended local binary patterns)
 729 //------------------------------------------------------------------------------
 730
 731 static Mat elbp(InputArray src, int radius, int neighbors) {
 732     Mat dst;
 733     elbp(src, dst, radius, neighbors);
 734     return dst;
 735 }
 736
 737 void LBPH::load(const FileStorage& fs) {
 738     fs["radius"] >> _radius;
 739     fs["neighbors"] >> _neighbors;
 740     fs["grid_x"] >> _grid_x;
 741     fs["grid_y"] >> _grid_y;
 742     //read matrices
 743     readFileNodeList(fs["histograms"], _histograms);
 744     fs["labels"] >> _labels;
 745 }
 746
 747 // See FaceRecognizer::save.
 748 void LBPH::save(FileStorage& fs) const {
 749     fs << "radius" << _radius;
 750     fs << "neighbors" << _neighbors;
 751     fs << "grid_x" << _grid_x;
 752     fs << "grid_y" << _grid_y;
 753     // write matrices
 754     writeFileNodeList(fs, "histograms", _histograms);
 755     fs << "labels" << _labels;
 756 }
 757
 758 void LBPH::train(InputArrayOfArrays _in_src, InputArray _in_labels) {
 759     this->train(_in_src, _in_labels, false);
 760 }
 761
 762 void LBPH::update(InputArrayOfArrays _in_src, InputArray _in_labels) {
 763     // got no data, just return
 764     if(_in_src.total() == 0)
 765         return;
 766
 767     this->train(_in_src, _in_labels, true);
 768 }
 769
 770 void LBPH::train(InputArrayOfArrays _in_src, InputArray _in_labels, bool preserveData) {
 771     if(_in_src.kind() != _InputArray::STD_VECTOR_MAT && _in_src.kind() != _InputArray::STD_VECTOR_VECTOR) {
 772         String error_message = "The images are expected as InputArray::STD_VECTOR_MAT (a std::vector<Mat>) or _InputArray::STD_VECTOR_VECTOR (a std::vector< std::vector<...> >).";
 773         CV_Error(Error::StsBadArg, error_message);
 774     }
 775     if(_in_src.total() == 0) {
 776         String error_message = format("Empty training data was given. You'll need more than one sample to learn a model.");
 777         CV_Error(Error::StsUnsupportedFormat, error_message);
 778     } else if(_in_labels.getMat().type() != CV_32SC1) {
 779         String error_message = format("Labels must be given as integer (CV_32SC1). Expected %d, but was %d.", CV_32SC1, _in_labels.type());
 780         CV_Error(Error::StsUnsupportedFormat, error_message);
 781     }
 782     // get the vector of matrices
 783     std::vector<Mat> src;
 784     _in_src.getMatVector(src);
 785     // get the label matrix
 786     Mat labels = _in_labels.getMat();
 787     // check if data is well- aligned
 788     if(labels.total() != src.size()) {
 789         String error_message = format("The number of samples (src) must equal the number of labels (labels). Was len(samples)=%d, len(labels)=%d.", src.size(), _labels.total());
 790         CV_Error(Error::StsBadArg, error_message);
 791     }
 792     // if this model should be trained without preserving old data, delete old model data
 793     if(!preserveData) {
 794         _labels.release();
 795         _histograms.clear();
 796     }
 797     // append labels to _labels matrix
 798     for(size_t labelIdx = 0; labelIdx < labels.total(); labelIdx++) {
 799         _labels.push_back(labels.at<int>((int)labelIdx));
 800     }
 801     // store the spatial histograms of the original data
 802     for(size_t sampleIdx = 0; sampleIdx < src.size(); sampleIdx++) {
 803         // calculate lbp image
 804         Mat lbp_image = elbp(src[sampleIdx], _radius, _neighbors);
 805         // get spatial histogram from this lbp image
 806         Mat p = spatial_histogram(
 807                 lbp_image, /* lbp_image */
 808                 static_cast<int>(std::pow(2.0, static_cast<double>(_neighbors))), /* number of possible patterns */
 809                 _grid_x, /* grid size x */
 810                 _grid_y, /* grid size y */
 811                 true);
 812         // add to templates
 813         _histograms.push_back(p);
 814     }
 815 }
 816
 817 void LBPH::predict(InputArray _src, int &minClass, double &minDist) const {
 818     if(_histograms.empty()) {
 819         // throw error if no data (or simply return -1?)
 820         String error_message = "This LBPH model is not computed yet. Did you call the train method?";
 821         CV_Error(Error::StsBadArg, error_message);
 822     }
 823     Mat src = _src.getMat();
 824     // get the spatial histogram from input image
 825     Mat lbp_image = elbp(src, _radius, _neighbors);
 826     Mat query = spatial_histogram(
 827             lbp_image, /* lbp_image */
 828             static_cast<int>(std::pow(2.0, static_cast<double>(_neighbors))), /* number of possible patterns */
 829             _grid_x, /* grid size x */
 830             _grid_y, /* grid size y */
 831             true /* normed histograms */);
 832     // find 1-nearest neighbor
 833     minDist = DBL_MAX;
 834     minClass = -1;
 835     for(size_t sampleIdx = 0; sampleIdx < _histograms.size(); sampleIdx++) {
 836         double dist = compareHist(_histograms[sampleIdx], query, HISTCMP_CHISQR_ALT);
 837         if((dist < minDist) && (dist < _threshold)) {
 838             minDist = dist;
 839             minClass = _labels.at<int>((int) sampleIdx);
 840         }
 841     }
 842 }
 843
 844 int LBPH::predict(InputArray _src) const {
 845     int label;
 846     double dummy;
 847     predict(_src, label, dummy);
 848     return label;
 849 }
 850
 851
 852 Ptr<FaceRecognizer> createEigenFaceRecognizer(int num_components, double threshold)
 853 {
 854     return makePtr<Eigenfaces>(num_components, threshold);
 855 }
 856
 857 Ptr<FaceRecognizer> createFisherFaceRecognizer(int num_components, double threshold)
 858 {
 859     return makePtr<Fisherfaces>(num_components, threshold);
 860 }
 861
 862 Ptr<FaceRecognizer> createLBPHFaceRecognizer(int radius, int neighbors,
 863                                              int grid_x, int grid_y, double threshold)
 864 {
 865     return makePtr<LBPH>(radius, neighbors, grid_x, grid_y, threshold);
 866 }
 867
 868 CV_INIT_ALGORITHM(Eigenfaces, "FaceRecognizer.Eigenfaces",
 869                   obj.info()->addParam(obj, "ncomponents", obj._num_components);
 870                   obj.info()->addParam(obj, "threshold", obj._threshold);
 871                   obj.info()->addParam(obj, "projections", obj._projections, true);
 872                   obj.info()->addParam(obj, "labels", obj._labels, true);
 873                   obj.info()->addParam(obj, "eigenvectors", obj._eigenvectors, true);
 874                   obj.info()->addParam(obj, "eigenvalues", obj._eigenvalues, true);
 875                   obj.info()->addParam(obj, "mean", obj._mean, true))
 876
 877 CV_INIT_ALGORITHM(Fisherfaces, "FaceRecognizer.Fisherfaces",
 878                   obj.info()->addParam(obj, "ncomponents", obj._num_components);
 879                   obj.info()->addParam(obj, "threshold", obj._threshold);
 880                   obj.info()->addParam(obj, "projections", obj._projections, true);
 881                   obj.info()->addParam(obj, "labels", obj._labels, true);
 882                   obj.info()->addParam(obj, "eigenvectors", obj._eigenvectors, true);
 883                   obj.info()->addParam(obj, "eigenvalues", obj._eigenvalues, true);
 884                   obj.info()->addParam(obj, "mean", obj._mean, true))
 885
 886 CV_INIT_ALGORITHM(LBPH, "FaceRecognizer.LBPH",
 887                   obj.info()->addParam(obj, "radius", obj._radius);
 888                   obj.info()->addParam(obj, "neighbors", obj._neighbors);
 889                   obj.info()->addParam(obj, "grid_x", obj._grid_x);
 890                   obj.info()->addParam(obj, "grid_y", obj._grid_y);
 891                   obj.info()->addParam(obj, "threshold", obj._threshold);
 892                   obj.info()->addParam(obj, "histograms", obj._histograms, true);
 893                   obj.info()->addParam(obj, "labels", obj._labels, true))
 894
 895 bool initModule_contrib()
 896 {
 897     Ptr<Algorithm> efaces = createEigenfaces_ptr_hidden(), ffaces = createFisherfaces_ptr_hidden(), lbph = createLBPH_ptr_hidden();
 898     return efaces->info() != 0 && ffaces->info() != 0 && lbph->info() != 0;
 899 }
 900
 901 }