opencv/src/cvaux/cvcalonder.cpp

   1 //*M///////////////////////////////////////////////////////////////////////////////////////
   2 //
   3 //  IMPORTANT: READ BEFORE DOWNLOADING, COPYING, INSTALLING OR USING.
   4 //
   5 //  By downloading, copying, installing or using the software you agree to this license.
   6 //  If you do not agree to this license, do not download, install,
   7 //  copy or use the software.
   8 //
   9 //
  10 //                           License Agreement
  11 //                For Open Source Computer Vision Library
  12 //
  13 // Copyright (C) 2000-2008, Intel Corporation, all rights reserved.
  14 // Copyright (C) 2009, Willow Garage Inc., all rights reserved.
  15 // Third party copyrights are property of their respective owners.
  16 //
  17 // Redistribution and use in source and binary forms, with or without modification,
  18 // are permitted provided that the following conditions are met:
  19 //
  20 //   * Redistribution's of source code must retain the above copyright notice,
  21 //     this list of conditions and the following disclaimer.
  22 //
  23 //   * Redistribution's in binary form must reproduce the above copyright notice,
  24 //     this list of conditions and the following disclaimer in the documentation
  25 //     and/or other materials provided with the distribution.
  26 //
  27 //   * The name of the copyright holders may not be used to endorse or promote products
  28 //     derived from this software without specific prior written permission.
  29 //
  30 // This software is provided by the copyright holders and contributors "as is" and
  31 // any express or implied warranties, including, but not limited to, the implied
  32 // warranties of merchantability and fitness for a particular purpose are disclaimed.
  33 // In no event shall the Intel Corporation or contributors be liable for any direct,
  34 // indirect, incidental, special, exemplary, or consequential damages
  35 // (including, but not limited to, procurement of substitute goods or services;
  36 // loss of use, data, or profits; or business interruption) however caused
  37 // and on any theory of liability, whether in contract, strict liability,
  38 // or tort (including negligence or otherwise) arising in any way out of
  39 // the use of this software, even if advised of the possibility of such damage.
  40 //
  41 //M*/
  42 #include "_cvaux.h"
  43
  44 #include <cxcore.h>
  45 #include <cvwimage.h>
  46 #include <vector>
  47 #include <iostream>
  48 #include <cv.h>
  49 #include <cmath>
  50 #include <cassert>
  51 #include <fstream>
  52 #include <cstring>
  53
  54 using namespace cv;
  55
  56
  57
  58 /****************************************************************************************\
  59 The code below is implementation of Calonder Descriptor and RTree Classifier
  60 originally introduced by Michael Calonder.
  61
  62 The code was integrated into OpenCV by Alexey Latyshev
  63 \****************************************************************************************/
  64
  65 namespace cv {
  66
  67
  68         //----------------------------
  69         //randomized_tree.cpp
  70
  71         inline uchar* getData(IplImage* image)
  72         {
  73                 return reinterpret_cast<uchar*>(image->imageData);
  74         }
  75
  76         inline float* RandomizedTree::getPosteriorByIndex(int index)
  77         {
  78                 return const_cast<float*>(const_cast<const RandomizedTree*>(this)->getPosteriorByIndex(index));
  79         }
  80
  81         inline const float* RandomizedTree::getPosteriorByIndex(int index) const
  82         {
  83                 return posteriors_[index];
  84         }
  85
  86         inline uchar* RandomizedTree::getPosteriorByIndex2(int index)
  87         {
  88                 return posteriors2_[index];
  89         }
  90
  91
  92         template < typename PointT >
  93         cv::WImageView1_b extractPatch(cv::WImageView1_b const& image, PointT pt, int patch_sz = PATCH_SIZE)
  94         {
  95                 const int offset = patch_sz / 2;
  96
  97                 // TODO: WImage{C}.View really should have const version
  98                 cv::WImageView1_b &img_ref = const_cast< cv::WImageView1_b& >(image);
  99                 return img_ref.View(pt.x - offset, pt.y - offset, patch_sz, patch_sz);
 100         }
 101
 102         template < typename PointT >
 103         cv::WImageView3_b extractPatch3(cv::WImageView3_b const& image, PointT pt)
 104         {
 105                 static const int offset = PATCH_SIZE / 2;
 106
 107                 // TODO: WImage{C}.View really should have const version
 108                 cv::WImageView3_b &img_ref = const_cast< cv::WImageView3_b& >(image);
 109                 return img_ref.View(pt.x - offset, pt.y - offset,
 110                         PATCH_SIZE, PATCH_SIZE);
 111         }
 112
 113         float *CSMatrixGenerator::cs_phi_   = NULL;
 114         int    CSMatrixGenerator::cs_phi_m_ = 0;
 115         int    CSMatrixGenerator::cs_phi_n_ = 0;
 116
 117         RandomizedTree::RandomizedTree()
 118                 : posteriors_(NULL), posteriors2_(NULL)
 119         {
 120         }
 121
 122         RandomizedTree::~RandomizedTree()
 123         {
 124                 freePosteriors(3);
 125         }
 126
 127         void RandomizedTree::createNodes(int num_nodes, cv::RNG &rng)
 128         {
 129                 nodes_.reserve(num_nodes);
 130                 for (int i = 0; i < num_nodes; ++i) {
 131                         nodes_.push_back( RTreeNode(rng(PATCH_SIZE),
 132                                 rng(PATCH_SIZE),
 133                                 rng(PATCH_SIZE),
 134                                 rng(PATCH_SIZE)) );
 135                 }
 136         }
 137
 138         int RandomizedTree::getIndex(uchar* patch_data) const
 139         {
 140                 int index = 0;
 141                 for (int d = 0; d < depth_; ++d) {
 142                         int child_offset = nodes_[index](patch_data);
 143                         index = 2*index + 1 + child_offset;
 144                 }
 145                 return index - nodes_.size();
 146         }
 147
 148         void RandomizedTree::train(std::vector<BaseKeypoint> const& base_set,
 149                 cv::RNG &rng, int depth, int views, size_t reduced_num_dim,
 150                 int num_quant_bits)
 151         {
 152
 153                 //CalonderPatchGenerator make_patch(NULL, rng);
 154                 PatchGenerator make_patch = PatchGenerator();
 155                 train(base_set, rng, make_patch, depth, views, reduced_num_dim, num_quant_bits);
 156         }
 157
 158         void RandomizedTree::train(std::vector<BaseKeypoint> const& base_set,
 159                 cv::RNG &rng, PatchGenerator &make_patch,
 160                 int depth, int views, size_t reduced_num_dim,
 161                 int num_quant_bits)
 162         {
 163                 init(base_set.size(), depth, rng);
 164
 165                 Mat patch;
 166
 167                 // Estimate posterior probabilities using random affine views
 168                 std::vector<BaseKeypoint>::const_iterator keypt_it;
 169                 int class_id = 0;
 170                 for (keypt_it = base_set.begin(); keypt_it != base_set.end(); ++keypt_it, ++class_id) {
 171                         for (int i = 0; i < views; ++i) {
 172
 173
 174                                 make_patch(keypt_it->image, Point2f(keypt_it->x,keypt_it->y) ,patch, Size(PATCH_SIZE,PATCH_SIZE),rng);
 175
 176                 IplImage _patch = patch;
 177                                 addExample(class_id, getData(&_patch));
 178                         }
 179                 }
 180
 181                 finalize(reduced_num_dim, num_quant_bits);
 182
 183         }
 184
 185         void RandomizedTree::allocPosteriorsAligned(int num_leaves, int num_classes)
 186         {
 187                 freePosteriors(3);
 188                 int err_cnt = 0;
 189
 190                 posteriors_ = new float*[num_leaves]; //(float**) malloc(num_leaves*sizeof(float*));
 191                 for (int i=0; i<num_leaves; ++i) {
 192                         //added
 193                         /* err_cnt += posix_memalign((void**)&posteriors_[i], 16, num_classes*sizeof(float));*/  posteriors_[i] = (float*)malloc(num_classes*sizeof(float));
 194                         memset(posteriors_[i], 0, num_classes*sizeof(float));
 195                 }
 196
 197                 posteriors2_ = new uchar*[num_leaves];
 198                 for (int i=0; i<num_leaves; ++i) {
 199                         //added
 200                         /*  err_cnt += posix_memalign((void**)&posteriors2_[i], 16, num_classes*sizeof(uchar)); */posteriors2_[i] = (uchar*)malloc(num_classes*sizeof(uchar));
 201                         memset(posteriors2_[i], 0, num_classes*sizeof(uchar));
 202                 }
 203
 204                 if (err_cnt) {
 205                         printf("Something went wrong in posix_memalign()! err_cnt=%i\n", err_cnt);
 206                         exit(0);
 207                 }
 208
 209                 classes_ = num_classes;
 210         }
 211
 212         void RandomizedTree::freePosteriors(int which)
 213         {
 214                 if (posteriors_ && (which&1)) {
 215                         for (int i=0; i<num_leaves_; ++i) {
 216                                 if (posteriors_[i]) {
 217                                         free(posteriors_[i]); //delete [] posteriors_[i];
 218                                         posteriors_[i] = NULL;
 219                                 }
 220                         }
 221                         delete [] posteriors_;
 222                         posteriors_ = NULL;
 223                 }
 224
 225                 if (posteriors2_ && (which&2)) {
 226                         for (int i=0; i<num_leaves_; ++i)
 227                                 free(posteriors2_[i]);
 228                         delete [] posteriors2_;
 229                         posteriors2_ = NULL;
 230                 }
 231
 232                 classes_ = -1;
 233         }
 234
 235         void RandomizedTree::init(int num_classes, int depth, cv::RNG &rng)
 236         {
 237                 depth_ = depth;
 238                 num_leaves_ = 1 << depth;        // 2**d
 239                 int num_nodes = num_leaves_ - 1; // 2**d - 1
 240
 241                 // Initialize probabilities and counts to 0
 242                 allocPosteriorsAligned(num_leaves_, num_classes);      // will set classes_ correctly
 243                 for (int i = 0; i < num_leaves_; ++i)
 244                         memset((void*)posteriors_[i], 0, num_classes*sizeof(float));
 245                 leaf_counts_.resize(num_leaves_);
 246
 247                 for (int i = 0; i < num_leaves_; ++i)
 248                         memset((void*)posteriors2_[i], 0, num_classes*sizeof(uchar));
 249
 250                 createNodes(num_nodes, rng);
 251         }
 252
 253         void RandomizedTree::addExample(int class_id, uchar* patch_data)
 254         {
 255                 int index = getIndex(patch_data);
 256                 float* posterior = getPosteriorByIndex(index);
 257                 ++leaf_counts_[index];
 258                 ++posterior[class_id];
 259         }
 260
 261         // returns the p% percentile of data (length n vector)
 262         static float percentile(float *data, int n, float p)
 263         {
 264                 assert(n>0);
 265                 assert(p>=0 && p<=1);
 266                 std::vector<float> vec(data, data+n);
 267                 sort(vec.begin(), vec.end());
 268                 int ix = (int)(p*(n-1));
 269                 return vec[ix];
 270         }
 271
 272         void RandomizedTree::finalize(size_t reduced_num_dim, int num_quant_bits)
 273         {
 274                 // Normalize by number of patches to reach each leaf
 275                 for (int index = 0; index < num_leaves_; ++index) {
 276                         float* posterior = posteriors_[index];
 277                         assert(posterior != NULL);
 278                         int count = leaf_counts_[index];
 279                         if (count != 0) {
 280                                 float normalizer = 1.0f / count;
 281                                 for (int c = 0; c < classes_; ++c) {
 282                                         *posterior *= normalizer;
 283                                         ++posterior;
 284                                 }
 285                         }
 286                 }
 287                 leaf_counts_.clear();
 288
 289                 // apply compressive sensing
 290                 if ((int)reduced_num_dim != classes_)
 291                         compressLeaves(reduced_num_dim);
 292                 else {
 293                         static bool notified = false;
 294                         //if (!notified)
 295                         //      printf("\n[OK] NO compression to leaves applied, dim=%i\n", reduced_num_dim);
 296                         notified = true;
 297                 }
 298
 299                 // convert float-posteriors to char-posteriors (quantization step)
 300                 makePosteriors2(num_quant_bits);
 301         }
 302
 303         void RandomizedTree::compressLeaves(size_t reduced_num_dim)
 304         {
 305                 static bool warned = false;
 306                 if (!warned) {
 307                         printf("\n[OK] compressing leaves with phi %i x %i\n", (int)reduced_num_dim, classes_);
 308                         warned = true;
 309                 }
 310
 311                 static bool warned2 = false;
 312                 if ((int)reduced_num_dim == classes_) {
 313                         if (!warned2)
 314                                 printf("[WARNING] RandomizedTree::compressLeaves:  not compressing because reduced_dim == classes()\n");
 315                         warned2 = true;
 316                         return;
 317                 }
 318
 319                 // DO NOT FREE RETURNED POINTER
 320                 float *cs_phi = CSMatrixGenerator::getCSMatrix(reduced_num_dim, classes_, CSMatrixGenerator::PDT_BERNOULLI);
 321
 322                 float *cs_posteriors = new float[num_leaves_ * reduced_num_dim];         // temp, num_leaves_ x reduced_num_dim
 323
 324                 for (int i=0; i<num_leaves_; ++i)
 325                 {
 326                         //added (inside cycle)
 327                         //float *post = getPosteriorByIndex(i);
 328                         //   float *prod = &cs_posteriors[i*reduced_num_dim];
 329                         //   cblas_sgemv(CblasRowMajor, CblasNoTrans, reduced_num_dim, classes_, 1.f, cs_phi,
 330                         //               classes_, post, 1, 0.f, prod, 1);
 331                         float *post = getPosteriorByIndex(i);
 332                         //Matrix multiplication
 333                         for (int idx = 0; idx < (int)reduced_num_dim; idx++)
 334                         {
 335                                 cs_posteriors[i*reduced_num_dim+idx] = 0.0f;
 336                                 for (int col = 0; col < classes_; col++)
 337                                 {
 338                                         cs_posteriors[i*reduced_num_dim+idx] += cs_phi[idx*reduced_num_dim + col] * post[col];
 339                                 }
 340                         }
 341                 }
 342
 343                 // copy new posteriors
 344                 freePosteriors(3);
 345                 allocPosteriorsAligned(num_leaves_, reduced_num_dim);
 346                 for (int i=0; i<num_leaves_; ++i)
 347                         memcpy(posteriors_[i], &cs_posteriors[i*reduced_num_dim], reduced_num_dim*sizeof(float));
 348                 classes_ = reduced_num_dim;
 349
 350                 delete [] cs_posteriors;
 351         }
 352
 353         void RandomizedTree::makePosteriors2(int num_quant_bits)
 354         {
 355                 int N = (1<<num_quant_bits) - 1;
 356
 357                 float perc[2];
 358                 estimateQuantPercForPosteriors(perc);
 359
 360                 assert(posteriors_ != NULL);
 361                 for (int i=0; i<num_leaves_; ++i)
 362                         quantizeVector(posteriors_[i], classes_, N, perc, posteriors2_[i]);
 363
 364                 // printf("makePosteriors2 quantization bounds: %.3e, %.3e (num_leaves=%i, N=%i)\n",
 365                 //        perc[0], perc[1], num_leaves_, N);
 366         }
 367
 368         void RandomizedTree::estimateQuantPercForPosteriors(float perc[2])
 369         {
 370                 // _estimate_ percentiles for this tree
 371                 // TODO: do this more accurately
 372                 assert(posteriors_ != NULL);
 373                 perc[0] = perc[1] = .0f;
 374                 for (int i=0; i<num_leaves_; i++) {
 375                         perc[0] += percentile(posteriors_[i], classes_, LOWER_QUANT_PERC);
 376                         perc[1] += percentile(posteriors_[i], classes_, UPPER_QUANT_PERC);
 377                 }
 378                 perc[0] /= num_leaves_;
 379                 perc[1] /= num_leaves_;
 380         }
 381
 382
 383         float* RandomizedTree::getPosterior(uchar* patch_data)
 384         {
 385                 return const_cast<float*>(const_cast<const RandomizedTree*>(this)->getPosterior(patch_data));
 386         }
 387
 388         const float* RandomizedTree::getPosterior(uchar* patch_data) const
 389         {
 390                 return getPosteriorByIndex( getIndex(patch_data) );
 391         }
 392
 393         uchar* RandomizedTree::getPosterior2(uchar* patch_data)
 394         {
 395                 return getPosteriorByIndex2( getIndex(patch_data) );
 396         }
 397
 398         void RandomizedTree::quantizeVector(float *vec, int dim, int N, float bnds[2], int clamp_mode)
 399         {
 400                 float map_bnd[2] = {0.f,(float)N};          // bounds of quantized target interval we're mapping to
 401                 for (int k=0; k<dim; ++k, ++vec) {
 402                         *vec = float(int((*vec - bnds[0])/(bnds[1] - bnds[0])*(map_bnd[1] - map_bnd[0]) + map_bnd[0]));
 403                         // 0: clamp both, lower and upper values
 404                         if (clamp_mode == 0)      *vec = (*vec<map_bnd[0]) ? map_bnd[0] : ((*vec>map_bnd[1]) ? map_bnd[1] : *vec);
 405                         // 1: clamp lower values only
 406                         else if (clamp_mode == 1) *vec = (*vec<map_bnd[0]) ? map_bnd[0] : *vec;
 407                         // 2: clamp upper values only
 408                         else if (clamp_mode == 2) *vec = (*vec>map_bnd[1]) ? map_bnd[1] : *vec;
 409                         // 4: no clamping
 410                         else if (clamp_mode == 4) ; // yep, nothing
 411                         else {
 412                                 printf("clamp_mode == %i is not valid (%s:%i).\n", clamp_mode, __FILE__, __LINE__);
 413                                 exit(1);
 414                         }
 415                 }
 416
 417         }
 418
 419         void RandomizedTree::quantizeVector(float *vec, int dim, int N, float bnds[2], uchar *dst)
 420         {
 421                 int map_bnd[2] = {0, N};          // bounds of quantized target interval we're mapping to
 422                 int tmp;
 423                 for (int k=0; k<dim; ++k) {
 424                         tmp = int((*vec - bnds[0])/(bnds[1] - bnds[0])*(map_bnd[1] - map_bnd[0]) + map_bnd[0]);
 425                         *dst = (uchar)((tmp<0) ? 0 : ((tmp>N) ? N : tmp));
 426                         ++vec;
 427                         ++dst;
 428                 }
 429         }
 430
 431
 432         void RandomizedTree::read(const char* file_name, int num_quant_bits)
 433         {
 434                 std::ifstream file(file_name, std::ifstream::binary);
 435                 read(file, num_quant_bits);
 436                 file.close();
 437         }
 438
 439         void RandomizedTree::read(std::istream &is, int num_quant_bits)
 440         {
 441                 is.read((char*)(&classes_), sizeof(classes_));
 442                 is.read((char*)(&depth_), sizeof(depth_));
 443
 444                 num_leaves_ = 1 << depth_;
 445                 int num_nodes = num_leaves_ - 1;
 446
 447                 nodes_.resize(num_nodes);
 448                 is.read((char*)(&nodes_[0]), num_nodes * sizeof(nodes_[0]));
 449
 450                 //posteriors_.resize(classes_ * num_leaves_);
 451                 //freePosteriors(3);
 452                 //printf("[DEBUG] reading: %i leaves, %i classes\n", num_leaves_, classes_);
 453                 allocPosteriorsAligned(num_leaves_, classes_);
 454                 for (int i=0; i<num_leaves_; i++)
 455                         is.read((char*)posteriors_[i], classes_ * sizeof(*posteriors_[0]));
 456
 457                 // make char-posteriors from float-posteriors
 458                 makePosteriors2(num_quant_bits);
 459         }
 460
 461         void RandomizedTree::write(const char* file_name) const
 462         {
 463                 std::ofstream file(file_name, std::ofstream::binary);
 464                 write(file);
 465                 file.close();
 466         }
 467
 468         void RandomizedTree::write(std::ostream &os) const
 469         {
 470                 if (!posteriors_) {
 471                         printf("WARNING: Cannot write float posteriors (posteriors_ = NULL).\n");
 472                         return;
 473                 }
 474
 475                 os.write((char*)(&classes_), sizeof(classes_));
 476                 os.write((char*)(&depth_), sizeof(depth_));
 477
 478                 os.write((char*)(&nodes_[0]), nodes_.size() * sizeof(nodes_[0]));
 479                 for (int i=0; i<num_leaves_; i++) {
 480                         os.write((char*)posteriors_[i], classes_ * sizeof(*posteriors_[0]));
 481                 }
 482         }
 483
 484
 485         void RandomizedTree::savePosteriors(std::string url, bool append)
 486         {
 487                 std::ofstream file(url.c_str(), (append?std::ios::app:std::ios::out));
 488                 for (int i=0; i<num_leaves_; i++) {
 489                         float *post = posteriors_[i];
 490                         char buf[20];
 491                         for (int i=0; i<classes_; i++) {
 492                                 sprintf(buf, "%.10e", *post++);
 493                                 file << buf << ((i<classes_-1) ? " " : "");
 494                         }
 495                         file << std::endl;
 496                 }
 497                 file.close();
 498         }
 499
 500         void RandomizedTree::savePosteriors2(std::string url, bool append)
 501         {
 502                 std::ofstream file(url.c_str(), (append?std::ios::app:std::ios::out));
 503                 for (int i=0; i<num_leaves_; i++) {
 504                         uchar *post = posteriors2_[i];
 505                         for (int i=0; i<classes_; i++)
 506                                 file << int(*post++) << (i<classes_-1?" ":"");
 507                         file << std::endl;
 508                 }
 509                 file.close();
 510         }
 511
 512         float* CSMatrixGenerator::getCSMatrix(int m, int n, PHI_DISTR_TYPE dt)
 513         {
 514                 assert(m <= n);
 515
 516                 if (cs_phi_m_!=m || cs_phi_n_!=n || cs_phi_==NULL) {
 517                         if (cs_phi_) delete [] cs_phi_;
 518                         cs_phi_ = new float[m*n];
 519                 }
 520
 521 #if 0 // debug - load the random matrix from a file (for reproducability of results)
 522                 //assert(m == 176);
 523                 //assert(n == 500);
 524                 //const char *phi = "/u/calonder/temp/dim_red/kpca_phi.txt";
 525                 const char *phi = "/u/calonder/temp/dim_red/debug_phi.txt";
 526                 std::ifstream ifs(phi);
 527                 for (size_t i=0; i<m*n; ++i) {
 528                         if (!ifs.good()) {
 529                                 printf("[ERROR] RandomizedTree::makeRandomMeasMatrix: problem reading '%s'\n", phi);
 530                                 exit(0);
 531                         }
 532                         ifs >> cs_phi[i];
 533                 }
 534                 ifs.close();
 535
 536                 static bool warned=false;
 537                 if (!warned) {
 538                         printf("[NOTE] RT: reading %ix%i PHI matrix from '%s'...\n", m, n, phi);
 539                         warned=true;
 540                 }
 541
 542                 return;
 543 #endif
 544
 545                 float *cs_phi = cs_phi_;
 546
 547                 if (m == n) {
 548                         // special case - set to 0 for safety
 549                         memset(cs_phi, 0, m*n*sizeof(float));
 550                         printf("[WARNING] %s:%i: square CS matrix (-> no reduction)\n", __FILE__, __LINE__);
 551                 }
 552                 else {
 553                         cv::RNG rng(23);
 554
 555                         // par is distr param, cf 'Favorable JL Distributions' (Baraniuk et al, 2006)
 556                         if (dt == PDT_GAUSS) {
 557                                 float par = (float)(1./m);
 558                                 //modified
 559                                 cv::RNG _rng;
 560                                 for (int i=0; i<m*n; ++i)
 561                                 {
 562                                         *cs_phi++ = (float)_rng.gaussian((double)par);//sample_normal<float>(0., par);
 563                                 }
 564                         }
 565                         else if (dt == PDT_BERNOULLI) {
 566                                 float par = (float)(1./sqrt((float)m));
 567                                 for (int i=0; i<m*n; ++i)
 568                                         *cs_phi++ = (rng(2)==0 ? par : -par);
 569                         }
 570                         else if (dt == PDT_DBFRIENDLY) {
 571                                 float par = (float)sqrt(3./m);
 572                                 for (int i=0; i<m*n; ++i) {
 573                                         //added
 574                                         int _i = rng(6);
 575                                         *cs_phi++ = (_i==0 ? par : (_i==1 ? -par : 0.f));
 576                                 }
 577                         }
 578                         else
 579                                 throw("PHI_DISTR_TYPE not implemented.");
 580                 }
 581
 582                 return cs_phi_;
 583         }
 584
 585         CSMatrixGenerator::~CSMatrixGenerator()
 586         {
 587                 if (cs_phi_) delete [] cs_phi_;
 588                 cs_phi_ = NULL;
 589         }
 590
 591
 592         //} // namespace features
 593
 594         //----------------------------
 595         //rtree_classifier.cpp
 596         //namespace features {
 597
 598         // Returns 16-byte aligned signatures that can be passed to getSignature().
 599         // Release by calling free() - NOT delete!
 600         //
 601         // note: 1) for num_sig>1 all signatures will still be 16-byte aligned, as
 602         //          long as sig_len%16 == 0 holds.
 603         //       2) casting necessary, otherwise it breaks gcc's strict aliasing rules
 604         inline void RTreeClassifier::safeSignatureAlloc(uchar **sig, int num_sig, int sig_len)
 605         {
 606                 assert(sig_len == 176);
 607                 void *p_sig;
 608                 //added
 609                 // posix_memalign(&p_sig, 16, num_sig*sig_len*sizeof(uchar));
 610                 p_sig = malloc(num_sig*sig_len*sizeof(uchar));
 611                 *sig = reinterpret_cast<uchar*>(p_sig);
 612         }
 613
 614         inline uchar* RTreeClassifier::safeSignatureAlloc(int num_sig, int sig_len)
 615         {
 616                 uchar *sig;
 617                 safeSignatureAlloc(&sig, num_sig, sig_len);
 618                 return sig;
 619         }
 620
 621         inline void add(int size, const float* src1, const float* src2, float* dst)
 622         {
 623                 while(--size >= 0) {
 624                         *dst = *src1 + *src2;
 625                         ++dst; ++src1; ++src2;
 626                 }
 627         }
 628
 629         inline void add(int size, const ushort* src1, const uchar* src2, ushort* dst)
 630         {
 631                 while(--size >= 0) {
 632                         *dst = *src1 + *src2;
 633                         ++dst; ++src1; ++src2;
 634                 }
 635         }
 636
 637         RTreeClassifier::RTreeClassifier()
 638                 : classes_(0)
 639         {
 640                 posteriors_ = NULL;
 641         }
 642
 643         void RTreeClassifier::train(std::vector<BaseKeypoint> const& base_set,
 644                 cv::RNG &rng, int num_trees, int depth,
 645                 int views, size_t reduced_num_dim,
 646                 int num_quant_bits, bool print_status)
 647         {
 648                 PatchGenerator make_patch = PatchGenerator();
 649                 train(base_set, rng, make_patch, num_trees, depth, views, reduced_num_dim, num_quant_bits, print_status);
 650         }
 651
 652         // Single-threaded version of train(), with progress output
 653         void RTreeClassifier::train(std::vector<BaseKeypoint> const& base_set,
 654                 cv::RNG &rng, PatchGenerator &make_patch, int num_trees,
 655                 int depth, int views, size_t reduced_num_dim,
 656                 int num_quant_bits, bool print_status)
 657         {
 658                 if (reduced_num_dim > base_set.size()) {
 659                         if (print_status)
 660                         {
 661                                 printf("INVALID PARAMS in RTreeClassifier::train: reduced_num_dim{%i} > base_set.size(){%i}\n",
 662                                         (int)reduced_num_dim, (int)base_set.size());
 663                         }
 664                         return;
 665                 }
 666
 667                 num_quant_bits_ = num_quant_bits;
 668                 classes_ = reduced_num_dim; // base_set.size();
 669                 original_num_classes_ = base_set.size();
 670                 trees_.resize(num_trees);
 671                 if (print_status)
 672                 {
 673                         printf("[OK] Training trees: base size=%i, reduced size=%i\n", (int)base_set.size(), (int)reduced_num_dim);
 674                 }
 675
 676                 int count = 1;
 677                 if (print_status)
 678                 {
 679                         printf("[OK] Trained 0 / %i trees", num_trees);  fflush(stdout);
 680                 }
 681                 //added
 682                 //BOOST_FOREACH( RandomizedTree &tree, trees_ ) {
 683                 //tree.train(base_set, rng, make_patch, depth, views, reduced_num_dim, num_quant_bits_);
 684                 //printf("\r[OK] Trained %i / %i trees", count++, num_trees);
 685                 //fflush(stdout);
 686                 for (int i=0; i<(int)trees_.size(); i++)
 687                 {
 688                         trees_[i].train(base_set, rng, make_patch, depth, views, reduced_num_dim, num_quant_bits_);
 689                         if (print_status)
 690                         {
 691                                 printf("\r[OK] Trained %i / %i trees", count++, num_trees);
 692                                 fflush(stdout);
 693                         }
 694                 }
 695
 696                 if (print_status)
 697                 {
 698                         printf("\n");
 699                         countZeroElements();
 700                         printf("\n\n");
 701                 }
 702         }
 703
 704         void RTreeClassifier::getSignature(IplImage* patch, float *sig)
 705         {
 706                 // Need pointer to 32x32 patch data
 707                 uchar buffer[PATCH_SIZE * PATCH_SIZE];
 708                 uchar* patch_data;
 709                 if (patch->widthStep != PATCH_SIZE) {
 710                         //printf("[INFO] patch is padded, data will be copied (%i/%i).\n",
 711                         //       patch->widthStep, PATCH_SIZE);
 712                         uchar* data = getData(patch);
 713                         patch_data = buffer;
 714                         for (int i = 0; i < PATCH_SIZE; ++i) {
 715                                 memcpy((void*)patch_data, (void*)data, PATCH_SIZE);
 716                                 data += patch->widthStep;
 717                                 patch_data += PATCH_SIZE;
 718                         }
 719                         patch_data = buffer;
 720                 }
 721                 else {
 722                         patch_data = getData(patch);
 723                 }
 724
 725                 memset((void*)sig, 0, classes_ * sizeof(float));
 726                 std::vector<RandomizedTree>::iterator tree_it;
 727
 728                 // get posteriors
 729                 float **posteriors = new float*[trees_.size()];  // TODO: move alloc outside this func
 730                 float **pp = posteriors;
 731                 for (tree_it = trees_.begin(); tree_it != trees_.end(); ++tree_it, pp++) {
 732                         *pp = tree_it->getPosterior(patch_data);
 733                         assert(*pp != NULL);
 734                 }
 735
 736                 // sum them up
 737                 pp = posteriors;
 738                 for (tree_it = trees_.begin(); tree_it != trees_.end(); ++tree_it, pp++)
 739                         add(classes_, sig, *pp, sig);
 740
 741                 delete [] posteriors;
 742                 posteriors = NULL;
 743
 744                 // full quantization (experimental)
 745 #if 0
 746                 int n_max = 1<<8 - 1;
 747                 int sum_max = (1<<4 - 1)*trees_.size();
 748                 int shift = 0;
 749                 while ((sum_max>>shift) > n_max) shift++;
 750
 751                 for (int i = 0; i < classes_; ++i) {
 752                         sig[i] = int(sig[i] + .5) >> shift;
 753                         if (sig[i]>n_max) sig[i] = n_max;
 754                 }
 755
 756                 static bool warned = false;
 757                 if (!warned) {
 758                         printf("[WARNING] Using full quantization (RTreeClassifier::getSignature)! shift=%i\n", shift);
 759                         warned = true;
 760                 }
 761 #else
 762                 // TODO: get rid of this multiply (-> number of trees is known at train
 763                 // time, exploit it in RandomizedTree::finalize())
 764                 float normalizer = 1.0f / trees_.size();
 765                 for (int i = 0; i < classes_; ++i)
 766                         sig[i] *= normalizer;
 767 #endif
 768         }
 769
 770
 771         // sum up 50 byte vectors of length 176
 772         // assume 5 bits max for input vector values
 773         // final shift is 3 bits right
 774         //void sum_50c_176c(uchar **pp, uchar *sig)
 775         //{
 776
 777         //}
 778
 779         void RTreeClassifier::getSignature(IplImage* patch, uchar *sig)
 780         {
 781                 // Need pointer to 32x32 patch data
 782                 uchar buffer[PATCH_SIZE * PATCH_SIZE];
 783                 uchar* patch_data;
 784                 if (patch->widthStep != PATCH_SIZE) {
 785                         //printf("[INFO] patch is padded, data will be copied (%i/%i).\n",
 786                         //       patch->widthStep, PATCH_SIZE);
 787                         uchar* data = getData(patch);
 788                         patch_data = buffer;
 789                         for (int i = 0; i < PATCH_SIZE; ++i) {
 790                                 memcpy((void*)patch_data, (void*)data, PATCH_SIZE);
 791                                 data += patch->widthStep;
 792                                 patch_data += PATCH_SIZE;
 793                         }
 794                         patch_data = buffer;
 795                 } else {
 796                         patch_data = getData(patch);
 797                 }
 798
 799                 std::vector<RandomizedTree>::iterator tree_it;
 800
 801                 // get posteriors
 802                 if (posteriors_ == NULL)
 803                 {
 804                         posteriors_ = new uchar*[trees_.size()];
 805                         //aadded
 806                         //  posix_memalign((void **)&ptemp_, 16, classes_*sizeof(ushort));
 807                         ptemp_ = (ushort*)malloc(classes_*sizeof(ushort));
 808                 }
 809                 uchar **pp = posteriors_;
 810                 for (tree_it = trees_.begin(); tree_it != trees_.end(); ++tree_it, pp++)
 811                         *pp = tree_it->getPosterior2(patch_data);
 812                 pp = posteriors_;
 813
 814 #if 0    // SSE2 optimized code
 815                 sum_50t_176c(pp, sig, ptemp_);    // sum them up
 816 #else
 817                 static bool warned = false;
 818
 819                 memset((void*)sig, 0, classes_ * sizeof(sig[0]));
 820                 ushort *sig16 = new ushort[classes_];           // TODO: make member, no alloc here
 821                 memset((void*)sig16, 0, classes_ * sizeof(sig16[0]));
 822                 for (tree_it = trees_.begin(); tree_it != trees_.end(); ++tree_it, pp++)
 823                         add(classes_, sig16, *pp, sig16);
 824
 825                 // squeeze signatures into an uchar
 826                 const bool full_shifting = true;
 827                 int shift;
 828                 if (full_shifting) {
 829                         float num_add_bits_f = log((float)trees_.size())/log(2.f);     // # additional bits required due to summation
 830                         int num_add_bits = int(num_add_bits_f);
 831                         if (num_add_bits_f != float(num_add_bits)) ++num_add_bits;
 832                         shift = num_quant_bits_ + num_add_bits - 8*sizeof(uchar);
 833                         //shift = num_quant_bits_ + num_add_bits - 2;
 834                         //shift = 6;
 835                         if (shift>0)
 836                                 for (int i = 0; i < classes_; ++i)
 837                                         sig[i] = (sig16[i] >> shift);      // &3 cut off all but lowest 2 bits, 3(dec) = 11(bin)
 838
 839                         if (!warned)
 840                                 printf("[OK] RTC: quantizing by FULL RIGHT SHIFT, shift = %i\n", shift);
 841                 }
 842                 else {
 843                         printf("[ERROR] RTC: not implemented!\n");
 844                         exit(0);
 845                 }
 846
 847                 if (!warned)
 848                         printf("[WARNING] RTC: unoptimized signature computation\n");
 849                 warned = true;
 850 #endif
 851         }
 852
 853
 854         void RTreeClassifier::getSparseSignature(IplImage *patch, float *sig, float thresh)
 855         {
 856                 getFloatSignature(patch, sig);
 857                 for (int i=0; i<classes_; ++i, sig++)
 858                         if (*sig < thresh) *sig = 0.f;
 859         }
 860
 861         int RTreeClassifier::countNonZeroElements(float *vec, int n, double tol)
 862         {
 863                 int res = 0;
 864                 while (n-- > 0)
 865                         res += (fabs(*vec++) > tol);
 866                 return res;
 867         }
 868
 869         void RTreeClassifier::read(const char* file_name)
 870         {
 871                 std::ifstream file(file_name, std::ifstream::binary);
 872                 read(file);
 873                 file.close();
 874         }
 875
 876         void RTreeClassifier::read(std::istream &is)
 877         {
 878                 int num_trees = 0;
 879                 is.read((char*)(&num_trees), sizeof(num_trees));
 880                 is.read((char*)(&classes_), sizeof(classes_));
 881                 is.read((char*)(&original_num_classes_), sizeof(original_num_classes_));
 882                 is.read((char*)(&num_quant_bits_), sizeof(num_quant_bits_));
 883
 884                 if (num_quant_bits_<1 || num_quant_bits_>8) {
 885                         printf("[WARNING] RTC: suspicious value num_quant_bits_=%i found; setting to %i.\n",
 886                                 num_quant_bits_, (int)DEFAULT_NUM_QUANT_BITS);
 887                         num_quant_bits_ = DEFAULT_NUM_QUANT_BITS;
 888                 }
 889
 890                 trees_.resize(num_trees);
 891                 std::vector<RandomizedTree>::iterator tree_it;
 892
 893                 for (tree_it = trees_.begin(); tree_it != trees_.end(); ++tree_it) {
 894                         tree_it->read(is, num_quant_bits_);
 895                 }
 896
 897                 printf("[OK] Loaded RTC, quantization=%i bits\n", num_quant_bits_);
 898
 899                 countZeroElements();
 900         }
 901
 902         void RTreeClassifier::write(const char* file_name) const
 903         {
 904                 std::ofstream file(file_name, std::ofstream::binary);
 905                 write(file);
 906                 file.close();
 907         }
 908
 909         void RTreeClassifier::write(std::ostream &os) const
 910         {
 911                 int num_trees = trees_.size();
 912                 os.write((char*)(&num_trees), sizeof(num_trees));
 913                 os.write((char*)(&classes_), sizeof(classes_));
 914                 os.write((char*)(&original_num_classes_), sizeof(original_num_classes_));
 915                 os.write((char*)(&num_quant_bits_), sizeof(num_quant_bits_));
 916                 printf("RTreeClassifier::write: num_quant_bits_=%i\n", num_quant_bits_);
 917
 918                 std::vector<RandomizedTree>::const_iterator tree_it;
 919                 for (tree_it = trees_.begin(); tree_it != trees_.end(); ++tree_it)
 920                         tree_it->write(os);
 921         }
 922
 923         void RTreeClassifier::saveAllFloatPosteriors(std::string url)
 924         {
 925                 printf("[DEBUG] writing all float posteriors to %s...\n", url.c_str());
 926                 for (int i=0; i<(int)trees_.size(); ++i)
 927                         trees_[i].savePosteriors(url, (i==0 ? false : true));
 928                 printf("[DEBUG] done\n");
 929         }
 930
 931         void RTreeClassifier::saveAllBytePosteriors(std::string url)
 932         {
 933                 printf("[DEBUG] writing all byte posteriors to %s...\n", url.c_str());
 934                 for (int i=0; i<(int)trees_.size(); ++i)
 935                         trees_[i].savePosteriors2(url, (i==0 ? false : true));
 936                 printf("[DEBUG] done\n");
 937         }
 938
 939
 940         void RTreeClassifier::setFloatPosteriorsFromTextfile_176(std::string url)
 941         {
 942                 std::ifstream ifs(url.c_str());
 943
 944                 for (int i=0; i<(int)trees_.size(); ++i) {
 945                         int num_classes = trees_[i].classes_;
 946                         assert(num_classes == 176);     // TODO: remove this limitation (arose due to SSE2 optimizations)
 947                         for (int k=0; k<trees_[i].num_leaves_; ++k) {
 948                                 float *post = trees_[i].getPosteriorByIndex(k);
 949                                 for (int j=0; j<num_classes; ++j, ++post)
 950                                         ifs >> *post;
 951                                 assert(ifs.good());
 952                         }
 953                 }
 954                 classes_ = 176;
 955
 956                 //setQuantization(num_quant_bits_);
 957
 958                 ifs.close();
 959                 printf("[EXPERIMENTAL] read entire tree from '%s'\n", url.c_str());
 960         }
 961
 962
 963         float RTreeClassifier::countZeroElements()
 964         {
 965                 int flt_zeros = 0;
 966                 int ui8_zeros = 0;
 967                 int num_elem = trees_[0].classes();
 968                 for (int i=0; i<(int)trees_.size(); ++i)
 969                         for (int k=0; k<(int)trees_[i].num_leaves_; ++k) {
 970                                 float *p = trees_[i].getPosteriorByIndex(k);
 971                                 uchar *p2 = trees_[i].getPosteriorByIndex2(k);
 972                                 assert(p); assert(p2);
 973                                 for (int j=0; j<num_elem; ++j, ++p, ++p2) {
 974                                         if (*p == 0.f) flt_zeros++;
 975                                         if (*p2 == 0) ui8_zeros++;
 976                                 }
 977                         }
 978                         num_elem = trees_.size()*trees_[0].num_leaves_*num_elem;
 979                         float flt_perc = 100.*flt_zeros/num_elem;
 980                         float ui8_perc = 100.*ui8_zeros/num_elem;
 981                         printf("[OK] RTC: overall %i/%i (%.3f%%) zeros in float leaves\n", flt_zeros, num_elem, flt_perc);
 982                         printf("          overall %i/%i (%.3f%%) zeros in uint8 leaves\n", ui8_zeros, num_elem, ui8_perc);
 983
 984                         return flt_perc;
 985         }
 986
 987         void RTreeClassifier::setQuantization(int num_quant_bits)
 988         {
 989                 for (int i=0; i<(int)trees_.size(); ++i)
 990                         trees_[i].applyQuantization(num_quant_bits);
 991
 992                 printf("[OK] signature quantization is now %i bits (before: %i)\n", num_quant_bits, num_quant_bits_);
 993                 num_quant_bits_ = num_quant_bits;
 994         }
 995
 996         void RTreeClassifier::discardFloatPosteriors()
 997         {
 998                 for (int i=0; i<(int)trees_.size(); ++i)
 999                         trees_[i].discardFloatPosteriors();
1000                 printf("[OK] RTC: discarded float posteriors of all trees\n");
1001         }
1002
1003         //} // namespace features
1004 }