apps/traincascade/boost.cpp

   1 #include "opencv2/core/core.hpp"
   2 #include "opencv2/core/internal.hpp"
   3
   4 #include "boost.h"
   5 #include "cascadeclassifier.h"
   6 #include <queue>
   7 #include "cxmisc.h"
   8
   9 using namespace std;
  10
  11 static inline double
  12 logRatio( double val )
  13 {
  14     const double eps = 1e-5;
  15
  16     val = max( val, eps );
  17     val = min( val, 1. - eps );
  18     return log( val/(1. - val) );
  19 }
  20
  21 #define CV_CMP_FLT(i,j) (i < j)
  22 static CV_IMPLEMENT_QSORT_EX( icvSortFlt, float, CV_CMP_FLT, const float* )
  23
  24 #define CV_CMP_NUM_IDX(i,j) (aux[i] < aux[j])
  25 static CV_IMPLEMENT_QSORT_EX( icvSortIntAux, int, CV_CMP_NUM_IDX, const float* )
  26 static CV_IMPLEMENT_QSORT_EX( icvSortUShAux, unsigned short, CV_CMP_NUM_IDX, const float* )
  27
  28 #define CV_THRESHOLD_EPS (0.00001F)
  29
  30 static const int MinBlockSize = 1 << 16;
  31 static const int BlockSizeDelta = 1 << 10;
  32
  33 // TODO remove this code duplication with ml/precomp.hpp
  34
  35 static int CV_CDECL icvCmpIntegers( const void* a, const void* b )
  36 {
  37     return *(const int*)a - *(const int*)b;
  38 }
  39
  40 static CvMat* cvPreprocessIndexArray( const CvMat* idx_arr, int data_arr_size, bool check_for_duplicates=false )
  41 {
  42     CvMat* idx = 0;
  43
  44     CV_FUNCNAME( "cvPreprocessIndexArray" );
  45
  46     __BEGIN__;
  47
  48     int i, idx_total, idx_selected = 0, step, type, prev = INT_MIN, is_sorted = 1;
  49     uchar* srcb = 0;
  50     int* srci = 0;
  51     int* dsti;
  52
  53     if( !CV_IS_MAT(idx_arr) )
  54         CV_ERROR( CV_StsBadArg, "Invalid index array" );
  55
  56     if( idx_arr->rows != 1 && idx_arr->cols != 1 )
  57         CV_ERROR( CV_StsBadSize, "the index array must be 1-dimensional" );
  58
  59     idx_total = idx_arr->rows + idx_arr->cols - 1;
  60     srcb = idx_arr->data.ptr;
  61     srci = idx_arr->data.i;
  62
  63     type = CV_MAT_TYPE(idx_arr->type);
  64     step = CV_IS_MAT_CONT(idx_arr->type) ? 1 : idx_arr->step/CV_ELEM_SIZE(type);
  65
  66     switch( type )
  67     {
  68     case CV_8UC1:
  69     case CV_8SC1:
  70         // idx_arr is array of 1's and 0's -
  71         // i.e. it is a mask of the selected components
  72         if( idx_total != data_arr_size )
  73             CV_ERROR( CV_StsUnmatchedSizes,
  74             "Component mask should contain as many elements as the total number of input variables" );
  75
  76         for( i = 0; i < idx_total; i++ )
  77             idx_selected += srcb[i*step] != 0;
  78
  79         if( idx_selected == 0 )
  80             CV_ERROR( CV_StsOutOfRange, "No components/input_variables is selected!" );
  81
  82         break;
  83     case CV_32SC1:
  84         // idx_arr is array of integer indices of selected components
  85         if( idx_total > data_arr_size )
  86             CV_ERROR( CV_StsOutOfRange,
  87             "index array may not contain more elements than the total number of input variables" );
  88         idx_selected = idx_total;
  89         // check if sorted already
  90         for( i = 0; i < idx_total; i++ )
  91         {
  92             int val = srci[i*step];
  93             if( val >= prev )
  94             {
  95                 is_sorted = 0;
  96                 break;
  97             }
  98             prev = val;
  99         }
 100         break;
 101     default:
 102         CV_ERROR( CV_StsUnsupportedFormat, "Unsupported index array data type "
 103                                            "(it should be 8uC1, 8sC1 or 32sC1)" );
 104     }
 105
 106     CV_CALL( idx = cvCreateMat( 1, idx_selected, CV_32SC1 ));
 107     dsti = idx->data.i;
 108
 109     if( type < CV_32SC1 )
 110     {
 111         for( i = 0; i < idx_total; i++ )
 112             if( srcb[i*step] )
 113                 *dsti++ = i;
 114     }
 115     else
 116     {
 117         for( i = 0; i < idx_total; i++ )
 118             dsti[i] = srci[i*step];
 119
 120         if( !is_sorted )
 121             qsort( dsti, idx_total, sizeof(dsti[0]), icvCmpIntegers );
 122
 123         if( dsti[0] < 0 || dsti[idx_total-1] >= data_arr_size )
 124             CV_ERROR( CV_StsOutOfRange, "the index array elements are out of range" );
 125
 126         if( check_for_duplicates )
 127         {
 128             for( i = 1; i < idx_total; i++ )
 129                 if( dsti[i] <= dsti[i-1] )
 130                     CV_ERROR( CV_StsBadArg, "There are duplicated index array elements" );
 131         }
 132     }
 133
 134     __END__;
 135
 136     if( cvGetErrStatus() < 0 )
 137         cvReleaseMat( &idx );
 138
 139     return idx;
 140 }
 141
 142 //----------------------------- CascadeBoostParams -------------------------------------------------
 143
 144 CvCascadeBoostParams::CvCascadeBoostParams() : minHitRate( 0.995F), maxFalseAlarm( 0.5F )
 145 {
 146     boost_type = CvBoost::GENTLE;
 147     use_surrogates = use_1se_rule = truncate_pruned_tree = false;
 148 }
 149
 150 CvCascadeBoostParams::CvCascadeBoostParams( int _boostType,
 151         float _minHitRate, float _maxFalseAlarm,
 152         double _weightTrimRate, int _maxDepth, int _maxWeakCount ) :
 153     CvBoostParams( _boostType, _maxWeakCount, _weightTrimRate, _maxDepth, false, 0 )
 154 {
 155     boost_type = CvBoost::GENTLE;
 156     minHitRate = _minHitRate;
 157     maxFalseAlarm = _maxFalseAlarm;
 158     use_surrogates = use_1se_rule = truncate_pruned_tree = false;
 159 }
 160
 161 void CvCascadeBoostParams::write( FileStorage &fs ) const
 162 {
 163     String boostTypeStr = boost_type == CvBoost::DISCRETE ? CC_DISCRETE_BOOST :
 164                           boost_type == CvBoost::REAL ? CC_REAL_BOOST :
 165                           boost_type == CvBoost::LOGIT ? CC_LOGIT_BOOST :
 166                           boost_type == CvBoost::GENTLE ? CC_GENTLE_BOOST : String();
 167     CV_Assert( !boostTypeStr.empty() );
 168     fs << CC_BOOST_TYPE << boostTypeStr;
 169     fs << CC_MINHITRATE << minHitRate;
 170     fs << CC_MAXFALSEALARM << maxFalseAlarm;
 171     fs << CC_TRIM_RATE << weight_trim_rate;
 172     fs << CC_MAX_DEPTH << max_depth;
 173     fs << CC_WEAK_COUNT << weak_count;
 174 }
 175
 176 bool CvCascadeBoostParams::read( const FileNode &node )
 177 {
 178     String boostTypeStr;
 179     FileNode rnode = node[CC_BOOST_TYPE];
 180     rnode >> boostTypeStr;
 181     boost_type = !boostTypeStr.compare( CC_DISCRETE_BOOST ) ? CvBoost::DISCRETE :
 182                  !boostTypeStr.compare( CC_REAL_BOOST ) ? CvBoost::REAL :
 183                  !boostTypeStr.compare( CC_LOGIT_BOOST ) ? CvBoost::LOGIT :
 184                  !boostTypeStr.compare( CC_GENTLE_BOOST ) ? CvBoost::GENTLE : -1;
 185     if (boost_type == -1)
 186         CV_Error( CV_StsBadArg, "unsupported Boost type" );
 187     node[CC_MINHITRATE] >> minHitRate;
 188     node[CC_MAXFALSEALARM] >> maxFalseAlarm;
 189     node[CC_TRIM_RATE] >> weight_trim_rate ;
 190     node[CC_MAX_DEPTH] >> max_depth ;
 191     node[CC_WEAK_COUNT] >> weak_count ;
 192     if ( minHitRate <= 0 || minHitRate > 1 ||
 193          maxFalseAlarm <= 0 || maxFalseAlarm > 1 ||
 194          weight_trim_rate <= 0 || weight_trim_rate > 1 ||
 195          max_depth <= 0 || weak_count <= 0 )
 196         CV_Error( CV_StsBadArg, "bad parameters range");
 197     return true;
 198 }
 199
 200 void CvCascadeBoostParams::printDefaults() const
 201 {
 202     cout << "--boostParams--" << endl;
 203     cout << "  [-bt <{" << CC_DISCRETE_BOOST << ", "
 204                         << CC_REAL_BOOST << ", "
 205                         << CC_LOGIT_BOOST ", "
 206                         << CC_GENTLE_BOOST << "(default)}>]" << endl;
 207     cout << "  [-minHitRate <min_hit_rate> = " << minHitRate << ">]" << endl;
 208     cout << "  [-maxFalseAlarmRate <max_false_alarm_rate = " << maxFalseAlarm << ">]" << endl;
 209     cout << "  [-weightTrimRate <weight_trim_rate = " << weight_trim_rate << ">]" << endl;
 210     cout << "  [-maxDepth <max_depth_of_weak_tree = " << max_depth << ">]" << endl;
 211     cout << "  [-maxWeakCount <max_weak_tree_count = " << weak_count << ">]" << endl;
 212 }
 213
 214 void CvCascadeBoostParams::printAttrs() const
 215 {
 216     String boostTypeStr = boost_type == CvBoost::DISCRETE ? CC_DISCRETE_BOOST :
 217                           boost_type == CvBoost::REAL ? CC_REAL_BOOST :
 218                           boost_type == CvBoost::LOGIT  ? CC_LOGIT_BOOST :
 219                           boost_type == CvBoost::GENTLE ? CC_GENTLE_BOOST : String();
 220     CV_Assert( !boostTypeStr.empty() );
 221     cout << "boostType: " << boostTypeStr << endl;
 222     cout << "minHitRate: " << minHitRate << endl;
 223     cout << "maxFalseAlarmRate: " <<  maxFalseAlarm << endl;
 224     cout << "weightTrimRate: " << weight_trim_rate << endl;
 225     cout << "maxDepth: " << max_depth << endl;
 226     cout << "maxWeakCount: " << weak_count << endl;
 227 }
 228
 229 bool CvCascadeBoostParams::scanAttr( const String prmName, const String val)
 230 {
 231     bool res = true;
 232
 233     if( !prmName.compare( "-bt" ) )
 234     {
 235         boost_type = !val.compare( CC_DISCRETE_BOOST ) ? CvBoost::DISCRETE :
 236                      !val.compare( CC_REAL_BOOST ) ? CvBoost::REAL :
 237                      !val.compare( CC_LOGIT_BOOST ) ? CvBoost::LOGIT :
 238                      !val.compare( CC_GENTLE_BOOST ) ? CvBoost::GENTLE : -1;
 239         if (boost_type == -1)
 240             res = false;
 241     }
 242     else if( !prmName.compare( "-minHitRate" ) )
 243     {
 244         minHitRate = (float) atof( val.c_str() );
 245     }
 246     else if( !prmName.compare( "-maxFalseAlarmRate" ) )
 247     {
 248         maxFalseAlarm = (float) atof( val.c_str() );
 249     }
 250     else if( !prmName.compare( "-weightTrimRate" ) )
 251     {
 252         weight_trim_rate = (float) atof( val.c_str() );
 253     }
 254     else if( !prmName.compare( "-maxDepth" ) )
 255     {
 256         max_depth = atoi( val.c_str() );
 257     }
 258     else if( !prmName.compare( "-maxWeakCount" ) )
 259     {
 260         weak_count = atoi( val.c_str() );
 261     }
 262     else
 263         res = false;
 264
 265     return res;
 266 }
 267
 268 CvDTreeNode* CvCascadeBoostTrainData::subsample_data( const CvMat* _subsample_idx )
 269 {
 270     CvDTreeNode* root = 0;
 271     CvMat* isubsample_idx = 0;
 272     CvMat* subsample_co = 0;
 273
 274     bool isMakeRootCopy = true;
 275
 276     if( !data_root )
 277         CV_Error( CV_StsError, "No training data has been set" );
 278
 279     if( _subsample_idx )
 280     {
 281         CV_Assert( (isubsample_idx = cvPreprocessIndexArray( _subsample_idx, sample_count )) != 0 );
 282
 283         if( isubsample_idx->cols + isubsample_idx->rows - 1 == sample_count )
 284         {
 285             const int* sidx = isubsample_idx->data.i;
 286             for( int i = 0; i < sample_count; i++ )
 287             {
 288                 if( sidx[i] != i )
 289                 {
 290                     isMakeRootCopy = false;
 291                     break;
 292                 }
 293             }
 294         }
 295         else
 296             isMakeRootCopy = false;
 297     }
 298
 299     if( isMakeRootCopy )
 300     {
 301         // make a copy of the root node
 302         CvDTreeNode temp;
 303         int i;
 304         root = new_node( 0, 1, 0, 0 );
 305         temp = *root;
 306         *root = *data_root;
 307         root->num_valid = temp.num_valid;
 308         if( root->num_valid )
 309         {
 310             for( i = 0; i < var_count; i++ )
 311                 root->num_valid[i] = data_root->num_valid[i];
 312         }
 313         root->cv_Tn = temp.cv_Tn;
 314         root->cv_node_risk = temp.cv_node_risk;
 315         root->cv_node_error = temp.cv_node_error;
 316     }
 317     else
 318     {
 319         int* sidx = isubsample_idx->data.i;
 320         // co - array of count/offset pairs (to handle duplicated values in _subsample_idx)
 321         int* co, cur_ofs = 0;
 322         int workVarCount = get_work_var_count();
 323         int count = isubsample_idx->rows + isubsample_idx->cols - 1;
 324
 325         root = new_node( 0, count, 1, 0 );
 326
 327         CV_Assert( (subsample_co = cvCreateMat( 1, sample_count*2, CV_32SC1 )) != 0);
 328         cvZero( subsample_co );
 329         co = subsample_co->data.i;
 330         for( int i = 0; i < count; i++ )
 331             co[sidx[i]*2]++;
 332         for( int i = 0; i < sample_count; i++ )
 333         {
 334             if( co[i*2] )
 335             {
 336                 co[i*2+1] = cur_ofs;
 337                 cur_ofs += co[i*2];
 338             }
 339             else
 340                 co[i*2+1] = -1;
 341         }
 342
 343         cv::AutoBuffer<uchar> inn_buf(sample_count*(2*sizeof(int) + sizeof(float)));
 344         // subsample ordered variables
 345         for( int vi = 0; vi < numPrecalcIdx; vi++ )
 346         {
 347             int ci = get_var_type(vi);
 348             CV_Assert( ci < 0 );
 349
 350             int *src_idx_buf = (int*)(uchar*)inn_buf;
 351             float *src_val_buf = (float*)(src_idx_buf + sample_count);
 352             int* sample_indices_buf = (int*)(src_val_buf + sample_count);
 353             const int* src_idx = 0;
 354             const float* src_val = 0;
 355             get_ord_var_data( data_root, vi, src_val_buf, src_idx_buf, &src_val, &src_idx, sample_indices_buf );
 356
 357             int j = 0, idx, count_i;
 358             int num_valid = data_root->get_num_valid(vi);
 359             CV_Assert( num_valid == sample_count );
 360
 361             if (is_buf_16u)
 362             {
 363                 unsigned short* udst_idx = (unsigned short*)(buf->data.s + root->buf_idx*get_length_subbuf() +
 364                     vi*sample_count + data_root->offset);
 365                 for( int i = 0; i < num_valid; i++ )
 366                 {
 367                     idx = src_idx[i];
 368                     count_i = co[idx*2];
 369                     if( count_i )
 370                         for( cur_ofs = co[idx*2+1]; count_i > 0; count_i--, j++, cur_ofs++ )
 371                             udst_idx[j] = (unsigned short)cur_ofs;
 372                 }
 373             }
 374             else
 375             {
 376                 int* idst_idx = buf->data.i + root->buf_idx*get_length_subbuf() +
 377                     vi*sample_count + root->offset;
 378                 for( int i = 0; i < num_valid; i++ )
 379                 {
 380                     idx = src_idx[i];
 381                     count_i = co[idx*2];
 382                     if( count_i )
 383                         for( cur_ofs = co[idx*2+1]; count_i > 0; count_i--, j++, cur_ofs++ )
 384                             idst_idx[j] = cur_ofs;
 385                 }
 386             }
 387         }
 388
 389         // subsample cv_lables
 390         const int* src_lbls = get_cv_labels(data_root, (int*)(uchar*)inn_buf);
 391         if (is_buf_16u)
 392         {
 393             unsigned short* udst = (unsigned short*)(buf->data.s + root->buf_idx*get_length_subbuf() +
 394                 (workVarCount-1)*sample_count + root->offset);
 395             for( int i = 0; i < count; i++ )
 396                 udst[i] = (unsigned short)src_lbls[sidx[i]];
 397         }
 398         else
 399         {
 400             int* idst = buf->data.i + root->buf_idx*get_length_subbuf() +
 401                 (workVarCount-1)*sample_count + root->offset;
 402             for( int i = 0; i < count; i++ )
 403                 idst[i] = src_lbls[sidx[i]];
 404         }
 405
 406         // subsample sample_indices
 407         const int* sample_idx_src = get_sample_indices(data_root, (int*)(uchar*)inn_buf);
 408         if (is_buf_16u)
 409         {
 410             unsigned short* sample_idx_dst = (unsigned short*)(buf->data.s + root->buf_idx*get_length_subbuf() +
 411                 workVarCount*sample_count + root->offset);
 412             for( int i = 0; i < count; i++ )
 413                 sample_idx_dst[i] = (unsigned short)sample_idx_src[sidx[i]];
 414         }
 415         else
 416         {
 417             int* sample_idx_dst = buf->data.i + root->buf_idx*get_length_subbuf() +
 418                 workVarCount*sample_count + root->offset;
 419             for( int i = 0; i < count; i++ )
 420                 sample_idx_dst[i] = sample_idx_src[sidx[i]];
 421         }
 422
 423         for( int vi = 0; vi < var_count; vi++ )
 424             root->set_num_valid(vi, count);
 425     }
 426
 427     cvReleaseMat( &isubsample_idx );
 428     cvReleaseMat( &subsample_co );
 429
 430     return root;
 431 }
 432
 433 //---------------------------- CascadeBoostTrainData -----------------------------
 434
 435 CvCascadeBoostTrainData::CvCascadeBoostTrainData( const CvFeatureEvaluator* _featureEvaluator,
 436                                                   const CvDTreeParams& _params )
 437 {
 438     is_classifier = true;
 439     var_all = var_count = (int)_featureEvaluator->getNumFeatures();
 440
 441     featureEvaluator = _featureEvaluator;
 442     shared = true;
 443     set_params( _params );
 444     max_c_count = MAX( 2, featureEvaluator->getMaxCatCount() );
 445     var_type = cvCreateMat( 1, var_count + 2, CV_32SC1 );
 446     if ( featureEvaluator->getMaxCatCount() > 0 )
 447     {
 448         numPrecalcIdx = 0;
 449         cat_var_count = var_count;
 450         ord_var_count = 0;
 451         for( int vi = 0; vi < var_count; vi++ )
 452         {
 453             var_type->data.i[vi] = vi;
 454         }
 455     }
 456     else
 457     {
 458         cat_var_count = 0;
 459         ord_var_count = var_count;
 460         for( int vi = 1; vi <= var_count; vi++ )
 461         {
 462             var_type->data.i[vi-1] = -vi;
 463         }
 464     }
 465     var_type->data.i[var_count] = cat_var_count;
 466     var_type->data.i[var_count+1] = cat_var_count+1;
 467
 468     int maxSplitSize = cvAlign(sizeof(CvDTreeSplit) + (MAX(0,max_c_count - 33)/32)*sizeof(int),sizeof(void*));
 469     int treeBlockSize = MAX((int)sizeof(CvDTreeNode)*8, maxSplitSize);
 470     treeBlockSize = MAX(treeBlockSize + BlockSizeDelta, MinBlockSize);
 471     tree_storage = cvCreateMemStorage( treeBlockSize );
 472     node_heap = cvCreateSet( 0, sizeof(node_heap[0]), sizeof(CvDTreeNode), tree_storage );
 473     split_heap = cvCreateSet( 0, sizeof(split_heap[0]), maxSplitSize, tree_storage );
 474 }
 475
 476 CvCascadeBoostTrainData::CvCascadeBoostTrainData( const CvFeatureEvaluator* _featureEvaluator,
 477                                                  int _numSamples,
 478                                                  int _precalcValBufSize, int _precalcIdxBufSize,
 479                                                  const CvDTreeParams& _params )
 480 {
 481     setData( _featureEvaluator, _numSamples, _precalcValBufSize, _precalcIdxBufSize, _params );
 482 }
 483
 484 void CvCascadeBoostTrainData::setData( const CvFeatureEvaluator* _featureEvaluator,
 485                                       int _numSamples,
 486                                       int _precalcValBufSize, int _precalcIdxBufSize,
 487                                       const CvDTreeParams& _params )
 488 {
 489     int* idst = 0;
 490     unsigned short* udst = 0;
 491
 492     uint64 effective_buf_size = 0;
 493     int effective_buf_height = 0, effective_buf_width = 0;
 494
 495
 496     clear();
 497     shared = true;
 498     have_labels = true;
 499     have_priors = false;
 500     is_classifier = true;
 501
 502     rng = &cv::theRNG();
 503
 504     set_params( _params );
 505
 506     CV_Assert( _featureEvaluator );
 507     featureEvaluator = _featureEvaluator;
 508
 509     max_c_count = MAX( 2, featureEvaluator->getMaxCatCount() );
 510     _resp = featureEvaluator->getCls();
 511     responses = &_resp;
 512     // TODO: check responses: elements must be 0 or 1
 513
 514     if( _precalcValBufSize < 0 || _precalcIdxBufSize < 0)
 515         CV_Error( CV_StsOutOfRange, "_numPrecalcVal and _numPrecalcIdx must be positive or 0" );
 516
 517     var_count = var_all = featureEvaluator->getNumFeatures() * featureEvaluator->getFeatureSize();
 518     sample_count = _numSamples;
 519
 520     is_buf_16u = false;
 521     if (sample_count < 65536)
 522         is_buf_16u = true;
 523
 524     numPrecalcVal = min( cvRound((double)_precalcValBufSize*1048576. / (sizeof(float)*sample_count)), var_count );
 525     numPrecalcIdx = min( cvRound((double)_precalcIdxBufSize*1048576. /
 526                 ((is_buf_16u ? sizeof(unsigned short) : sizeof (int))*sample_count)), var_count );
 527
 528     assert( numPrecalcIdx >= 0 && numPrecalcVal >= 0 );
 529
 530     valCache.create( numPrecalcVal, sample_count, CV_32FC1 );
 531     var_type = cvCreateMat( 1, var_count + 2, CV_32SC1 );
 532
 533     if ( featureEvaluator->getMaxCatCount() > 0 )
 534     {
 535         numPrecalcIdx = 0;
 536         cat_var_count = var_count;
 537         ord_var_count = 0;
 538         for( int vi = 0; vi < var_count; vi++ )
 539         {
 540             var_type->data.i[vi] = vi;
 541         }
 542     }
 543     else
 544     {
 545         cat_var_count = 0;
 546         ord_var_count = var_count;
 547         for( int vi = 1; vi <= var_count; vi++ )
 548         {
 549             var_type->data.i[vi-1] = -vi;
 550         }
 551     }
 552     var_type->data.i[var_count] = cat_var_count;
 553     var_type->data.i[var_count+1] = cat_var_count+1;
 554     work_var_count = ( cat_var_count ? 0 : numPrecalcIdx ) + 1/*cv_lables*/;
 555     buf_count = 2;
 556
 557     buf_size = -1; // the member buf_size is obsolete
 558
 559     effective_buf_size = (uint64)(work_var_count + 1)*(uint64)sample_count * buf_count; // this is the total size of "CvMat buf" to be allocated
 560     effective_buf_width = sample_count;
 561     effective_buf_height = work_var_count+1;
 562
 563     if (effective_buf_width >= effective_buf_height)
 564         effective_buf_height *= buf_count;
 565     else
 566         effective_buf_width *= buf_count;
 567
 568     if ((uint64)effective_buf_width * (uint64)effective_buf_height != effective_buf_size)
 569     {
 570         CV_Error(CV_StsBadArg, "The memory buffer cannot be allocated since its size exceeds integer fields limit");
 571     }
 572
 573     if ( is_buf_16u )
 574         buf = cvCreateMat( effective_buf_height, effective_buf_width, CV_16UC1 );
 575     else
 576         buf = cvCreateMat( effective_buf_height, effective_buf_width, CV_32SC1 );
 577
 578     cat_count = cvCreateMat( 1, cat_var_count + 1, CV_32SC1 );
 579
 580     // precalculate valCache and set indices in buf
 581     precalculate();
 582
 583     // now calculate the maximum size of split,
 584     // create memory storage that will keep nodes and splits of the decision tree
 585     // allocate root node and the buffer for the whole training data
 586     int maxSplitSize = cvAlign(sizeof(CvDTreeSplit) +
 587         (MAX(0,sample_count - 33)/32)*sizeof(int),sizeof(void*));
 588     int treeBlockSize = MAX((int)sizeof(CvDTreeNode)*8, maxSplitSize);
 589     treeBlockSize = MAX(treeBlockSize + BlockSizeDelta, MinBlockSize);
 590     tree_storage = cvCreateMemStorage( treeBlockSize );
 591     node_heap = cvCreateSet( 0, sizeof(*node_heap), sizeof(CvDTreeNode), tree_storage );
 592
 593     int nvSize = var_count*sizeof(int);
 594     nvSize = cvAlign(MAX( nvSize, (int)sizeof(CvSetElem) ), sizeof(void*));
 595     int tempBlockSize = nvSize;
 596     tempBlockSize = MAX( tempBlockSize + BlockSizeDelta, MinBlockSize );
 597     temp_storage = cvCreateMemStorage( tempBlockSize );
 598     nv_heap = cvCreateSet( 0, sizeof(*nv_heap), nvSize, temp_storage );
 599
 600     data_root = new_node( 0, sample_count, 0, 0 );
 601
 602     // set sample labels
 603     if (is_buf_16u)
 604         udst = (unsigned short*)(buf->data.s + work_var_count*sample_count);
 605     else
 606         idst = buf->data.i + work_var_count*sample_count;
 607
 608     for (int si = 0; si < sample_count; si++)
 609     {
 610         if (udst)
 611             udst[si] = (unsigned short)si;
 612         else
 613             idst[si] = si;
 614     }
 615     for( int vi = 0; vi < var_count; vi++ )
 616         data_root->set_num_valid(vi, sample_count);
 617     for( int vi = 0; vi < cat_var_count; vi++ )
 618         cat_count->data.i[vi] = max_c_count;
 619
 620     cat_count->data.i[cat_var_count] = 2;
 621
 622     maxSplitSize = cvAlign(sizeof(CvDTreeSplit) +
 623         (MAX(0,max_c_count - 33)/32)*sizeof(int),sizeof(void*));
 624     split_heap = cvCreateSet( 0, sizeof(*split_heap), maxSplitSize, tree_storage );
 625
 626     priors = cvCreateMat( 1, get_num_classes(), CV_64F );
 627     cvSet(priors, cvScalar(1));
 628     priors_mult = cvCloneMat( priors );
 629     counts = cvCreateMat( 1, get_num_classes(), CV_32SC1 );
 630     direction = cvCreateMat( 1, sample_count, CV_8UC1 );
 631     split_buf = cvCreateMat( 1, sample_count, CV_32SC1 );//TODO: make a pointer
 632 }
 633
 634 void CvCascadeBoostTrainData::free_train_data()
 635 {
 636     CvDTreeTrainData::free_train_data();
 637     valCache.release();
 638 }
 639
 640 const int* CvCascadeBoostTrainData::get_class_labels( CvDTreeNode* n, int* labelsBuf)
 641 {
 642     int nodeSampleCount = n->sample_count;
 643     int rStep = CV_IS_MAT_CONT( responses->type ) ? 1 : responses->step / CV_ELEM_SIZE( responses->type );
 644
 645     int* sampleIndicesBuf = labelsBuf; //
 646     const int* sampleIndices = get_sample_indices(n, sampleIndicesBuf);
 647     for( int si = 0; si < nodeSampleCount; si++ )
 648     {
 649         int sidx = sampleIndices[si];
 650         labelsBuf[si] = (int)responses->data.fl[sidx*rStep];
 651     }
 652     return labelsBuf;
 653 }
 654
 655 const int* CvCascadeBoostTrainData::get_sample_indices( CvDTreeNode* n, int* indicesBuf )
 656 {
 657     return CvDTreeTrainData::get_cat_var_data( n, get_work_var_count(), indicesBuf );
 658 }
 659
 660 const int* CvCascadeBoostTrainData::get_cv_labels( CvDTreeNode* n, int* labels_buf )
 661 {
 662     return CvDTreeTrainData::get_cat_var_data( n, get_work_var_count() - 1, labels_buf );
 663 }
 664
 665 void CvCascadeBoostTrainData::get_ord_var_data( CvDTreeNode* n, int vi, float* ordValuesBuf, int* sortedIndicesBuf,
 666         const float** ordValues, const int** sortedIndices, int* sampleIndicesBuf )
 667 {
 668     int nodeSampleCount = n->sample_count;
 669     const int* sampleIndices = get_sample_indices(n, sampleIndicesBuf);
 670
 671     if ( vi < numPrecalcIdx )
 672     {
 673         if( !is_buf_16u )
 674             *sortedIndices = buf->data.i + n->buf_idx*get_length_subbuf() + vi*sample_count + n->offset;
 675         else
 676         {
 677             const unsigned short* shortIndices = (const unsigned short*)(buf->data.s + n->buf_idx*get_length_subbuf() +
 678                                                     vi*sample_count + n->offset );
 679             for( int i = 0; i < nodeSampleCount; i++ )
 680                 sortedIndicesBuf[i] = shortIndices[i];
 681
 682             *sortedIndices = sortedIndicesBuf;
 683         }
 684
 685         if( vi < numPrecalcVal )
 686         {
 687             for( int i = 0; i < nodeSampleCount; i++ )
 688             {
 689                 int idx = (*sortedIndices)[i];
 690                 idx = sampleIndices[idx];
 691                 ordValuesBuf[i] =  valCache.at<float>( vi, idx);
 692             }
 693         }
 694         else
 695         {
 696             for( int i = 0; i < nodeSampleCount; i++ )
 697             {
 698                 int idx = (*sortedIndices)[i];
 699                 idx = sampleIndices[idx];
 700                 ordValuesBuf[i] = (*featureEvaluator)( vi, idx);
 701             }
 702         }
 703     }
 704     else // vi >= numPrecalcIdx
 705     {
 706         cv::AutoBuffer<float> abuf(nodeSampleCount);
 707         float* sampleValues = &abuf[0];
 708
 709         if ( vi < numPrecalcVal )
 710         {
 711             for( int i = 0; i < nodeSampleCount; i++ )
 712             {
 713                 sortedIndicesBuf[i] = i;
 714                 sampleValues[i] = valCache.at<float>( vi, sampleIndices[i] );
 715             }
 716         }
 717         else
 718         {
 719             for( int i = 0; i < nodeSampleCount; i++ )
 720             {
 721                 sortedIndicesBuf[i] = i;
 722                 sampleValues[i] = (*featureEvaluator)( vi, sampleIndices[i]);
 723             }
 724         }
 725         icvSortIntAux( sortedIndicesBuf, nodeSampleCount, &sampleValues[0] );
 726         for( int i = 0; i < nodeSampleCount; i++ )
 727             ordValuesBuf[i] = (&sampleValues[0])[sortedIndicesBuf[i]];
 728         *sortedIndices = sortedIndicesBuf;
 729     }
 730
 731     *ordValues = ordValuesBuf;
 732 }
 733
 734 const int* CvCascadeBoostTrainData::get_cat_var_data( CvDTreeNode* n, int vi, int* catValuesBuf )
 735 {
 736     int nodeSampleCount = n->sample_count;
 737     int* sampleIndicesBuf = catValuesBuf; //
 738     const int* sampleIndices = get_sample_indices(n, sampleIndicesBuf);
 739
 740     if ( vi < numPrecalcVal )
 741     {
 742         for( int i = 0; i < nodeSampleCount; i++ )
 743             catValuesBuf[i] = (int) valCache.at<float>( vi, sampleIndices[i]);
 744     }
 745     else
 746     {
 747         if( vi >= numPrecalcVal && vi < var_count )
 748         {
 749             for( int i = 0; i < nodeSampleCount; i++ )
 750                 catValuesBuf[i] = (int)(*featureEvaluator)( vi, sampleIndices[i] );
 751         }
 752         else
 753         {
 754             get_cv_labels( n, catValuesBuf );
 755         }
 756     }
 757
 758     return catValuesBuf;
 759 }
 760
 761 float CvCascadeBoostTrainData::getVarValue( int vi, int si )
 762 {
 763     if ( vi < numPrecalcVal && !valCache.empty() )
 764         return valCache.at<float>( vi, si );
 765     return (*featureEvaluator)( vi, si );
 766 }
 767
 768
 769 struct FeatureIdxOnlyPrecalc
 770 {
 771     FeatureIdxOnlyPrecalc( const CvFeatureEvaluator* _featureEvaluator, CvMat* _buf, int _sample_count, bool _is_buf_16u )
 772     {
 773         featureEvaluator = _featureEvaluator;
 774         sample_count = _sample_count;
 775         udst = (unsigned short*)_buf->data.s;
 776         idst = _buf->data.i;
 777         is_buf_16u = _is_buf_16u;
 778     }
 779     void operator()( const BlockedRange& range ) const
 780     {
 781         cv::AutoBuffer<float> valCache(sample_count);
 782         float* valCachePtr = (float*)valCache;
 783         for ( int fi = range.begin(); fi < range.end(); fi++)
 784         {
 785             for( int si = 0; si < sample_count; si++ )
 786             {
 787                 valCachePtr[si] = (*featureEvaluator)( fi, si );
 788                 if ( is_buf_16u )
 789                     *(udst + fi*sample_count + si) = (unsigned short)si;
 790                 else
 791                     *(idst + fi*sample_count + si) = si;
 792             }
 793             if ( is_buf_16u )
 794                 icvSortUShAux( udst + fi*sample_count, sample_count, valCachePtr );
 795             else
 796                 icvSortIntAux( idst + fi*sample_count, sample_count, valCachePtr );
 797         }
 798     }
 799     const CvFeatureEvaluator* featureEvaluator;
 800     int sample_count;
 801     int* idst;
 802     unsigned short* udst;
 803     bool is_buf_16u;
 804 };
 805
 806 struct FeatureValAndIdxPrecalc
 807 {
 808     FeatureValAndIdxPrecalc( const CvFeatureEvaluator* _featureEvaluator, CvMat* _buf, Mat* _valCache, int _sample_count, bool _is_buf_16u )
 809     {
 810         featureEvaluator = _featureEvaluator;
 811         valCache = _valCache;
 812         sample_count = _sample_count;
 813         udst = (unsigned short*)_buf->data.s;
 814         idst = _buf->data.i;
 815         is_buf_16u = _is_buf_16u;
 816     }
 817     void operator()( const BlockedRange& range ) const
 818     {
 819         for ( int fi = range.begin(); fi < range.end(); fi++)
 820         {
 821             for( int si = 0; si < sample_count; si++ )
 822             {
 823                 valCache->at<float>(fi,si) = (*featureEvaluator)( fi, si );
 824                 if ( is_buf_16u )
 825                     *(udst + fi*sample_count + si) = (unsigned short)si;
 826                 else
 827                     *(idst + fi*sample_count + si) = si;
 828             }
 829             if ( is_buf_16u )
 830                 icvSortUShAux( udst + fi*sample_count, sample_count, valCache->ptr<float>(fi) );
 831             else
 832                 icvSortIntAux( idst + fi*sample_count, sample_count, valCache->ptr<float>(fi) );
 833         }
 834     }
 835     const CvFeatureEvaluator* featureEvaluator;
 836     Mat* valCache;
 837     int sample_count;
 838     int* idst;
 839     unsigned short* udst;
 840     bool is_buf_16u;
 841 };
 842
 843 struct FeatureValOnlyPrecalc
 844 {
 845     FeatureValOnlyPrecalc( const CvFeatureEvaluator* _featureEvaluator, Mat* _valCache, int _sample_count )
 846     {
 847         featureEvaluator = _featureEvaluator;
 848         valCache = _valCache;
 849         sample_count = _sample_count;
 850     }
 851     void operator()( const BlockedRange& range ) const
 852     {
 853         for ( int fi = range.begin(); fi < range.end(); fi++)
 854             for( int si = 0; si < sample_count; si++ )
 855                 valCache->at<float>(fi,si) = (*featureEvaluator)( fi, si );
 856     }
 857     const CvFeatureEvaluator* featureEvaluator;
 858     Mat* valCache;
 859     int sample_count;
 860 };
 861
 862 void CvCascadeBoostTrainData::precalculate()
 863 {
 864     int minNum = MIN( numPrecalcVal, numPrecalcIdx);
 865
 866     double proctime = -TIME( 0 );
 867     parallel_for( BlockedRange(numPrecalcVal, numPrecalcIdx),
 868                   FeatureIdxOnlyPrecalc(featureEvaluator, buf, sample_count, is_buf_16u!=0) );
 869     parallel_for( BlockedRange(0, minNum),
 870                   FeatureValAndIdxPrecalc(featureEvaluator, buf, &valCache, sample_count, is_buf_16u!=0) );
 871     parallel_for( BlockedRange(minNum, numPrecalcVal),
 872                   FeatureValOnlyPrecalc(featureEvaluator, &valCache, sample_count) );
 873     cout << "Precalculation time: " << (proctime + TIME( 0 )) << endl;
 874 }
 875
 876 //-------------------------------- CascadeBoostTree ----------------------------------------
 877
 878 CvDTreeNode* CvCascadeBoostTree::predict( int sampleIdx ) const
 879 {
 880     CvDTreeNode* node = root;
 881     if( !node )
 882         CV_Error( CV_StsError, "The tree has not been trained yet" );
 883
 884     if ( ((CvCascadeBoostTrainData*)data)->featureEvaluator->getMaxCatCount() == 0 ) // ordered
 885     {
 886         while( node->left )
 887         {
 888             CvDTreeSplit* split = node->split;
 889             float val = ((CvCascadeBoostTrainData*)data)->getVarValue( split->var_idx, sampleIdx );
 890             node = val <= split->ord.c ? node->left : node->right;
 891         }
 892     }
 893     else // categorical
 894     {
 895         while( node->left )
 896         {
 897             CvDTreeSplit* split = node->split;
 898             int c = (int)((CvCascadeBoostTrainData*)data)->getVarValue( split->var_idx, sampleIdx );
 899             node = CV_DTREE_CAT_DIR(c, split->subset) < 0 ? node->left : node->right;
 900         }
 901     }
 902     return node;
 903 }
 904
 905 void CvCascadeBoostTree::write( FileStorage &fs, const Mat& featureMap )
 906 {
 907     int maxCatCount = ((CvCascadeBoostTrainData*)data)->featureEvaluator->getMaxCatCount();
 908     int subsetN = (maxCatCount + 31)/32;
 909     queue<CvDTreeNode*> internalNodesQueue;
 910     int size = (int)pow( 2.f, (float)ensemble->get_params().max_depth);
 911     Ptr<float> leafVals = new float[size];
 912     int leafValIdx = 0;
 913     int internalNodeIdx = 1;
 914     CvDTreeNode* tempNode;
 915
 916     CV_DbgAssert( root );
 917     internalNodesQueue.push( root );
 918
 919     fs << "{";
 920     fs << CC_INTERNAL_NODES << "[:";
 921     while (!internalNodesQueue.empty())
 922     {
 923         tempNode = internalNodesQueue.front();
 924         CV_Assert( tempNode->left );
 925         if ( !tempNode->left->left && !tempNode->left->right) // left node is leaf
 926         {
 927             leafVals[-leafValIdx] = (float)tempNode->left->value;
 928             fs << leafValIdx-- ;
 929         }
 930         else
 931         {
 932             internalNodesQueue.push( tempNode->left );
 933             fs << internalNodeIdx++;
 934         }
 935         CV_Assert( tempNode->right );
 936         if ( !tempNode->right->left && !tempNode->right->right) // right node is leaf
 937         {
 938             leafVals[-leafValIdx] = (float)tempNode->right->value;
 939             fs << leafValIdx--;
 940         }
 941         else
 942         {
 943             internalNodesQueue.push( tempNode->right );
 944             fs << internalNodeIdx++;
 945         }
 946         int fidx = tempNode->split->var_idx;
 947         fidx = featureMap.empty() ? fidx : featureMap.at<int>(0, fidx);
 948         fs << fidx;
 949         if ( !maxCatCount )
 950             fs << tempNode->split->ord.c;
 951         else
 952             for( int i = 0; i < subsetN; i++ )
 953                 fs << tempNode->split->subset[i];
 954         internalNodesQueue.pop();
 955     }
 956     fs << "]"; // CC_INTERNAL_NODES
 957
 958     fs << CC_LEAF_VALUES << "[:";
 959     for (int ni = 0; ni < -leafValIdx; ni++)
 960         fs << leafVals[ni];
 961     fs << "]"; // CC_LEAF_VALUES
 962     fs << "}";
 963 }
 964
 965 void CvCascadeBoostTree::read( const FileNode &node, CvBoost* _ensemble,
 966                                 CvDTreeTrainData* _data )
 967 {
 968     int maxCatCount = ((CvCascadeBoostTrainData*)_data)->featureEvaluator->getMaxCatCount();
 969     int subsetN = (maxCatCount + 31)/32;
 970     int step = 3 + ( maxCatCount>0 ? subsetN : 1 );
 971
 972     queue<CvDTreeNode*> internalNodesQueue;
 973     FileNodeIterator internalNodesIt, leafValsuesIt;
 974     CvDTreeNode* prntNode, *cldNode;
 975
 976     clear();
 977     data = _data;
 978     ensemble = _ensemble;
 979     pruned_tree_idx = 0;
 980
 981     // read tree nodes
 982     FileNode rnode = node[CC_INTERNAL_NODES];
 983     internalNodesIt = rnode.end();
 984     leafValsuesIt = node[CC_LEAF_VALUES].end();
 985     internalNodesIt--; leafValsuesIt--;
 986     for( size_t i = 0; i < rnode.size()/step; i++ )
 987     {
 988         prntNode = data->new_node( 0, 0, 0, 0 );
 989         if ( maxCatCount > 0 )
 990         {
 991             prntNode->split = data->new_split_cat( 0, 0 );
 992             for( int j = subsetN-1; j>=0; j--)
 993             {
 994                 *internalNodesIt >> prntNode->split->subset[j]; internalNodesIt--;
 995             }
 996         }
 997         else
 998         {
 999             float split_value;
1000             *internalNodesIt >> split_value; internalNodesIt--;
1001             prntNode->split = data->new_split_ord( 0, split_value, 0, 0, 0);
1002         }
1003         *internalNodesIt >> prntNode->split->var_idx; internalNodesIt--;
1004         int ridx, lidx;
1005         *internalNodesIt >> ridx; internalNodesIt--;
1006         *internalNodesIt >> lidx;internalNodesIt--;
1007         if ( ridx <= 0)
1008         {
1009             prntNode->right = cldNode = data->new_node( 0, 0, 0, 0 );
1010             *leafValsuesIt >> cldNode->value; leafValsuesIt--;
1011             cldNode->parent = prntNode;
1012         }
1013         else
1014         {
1015             prntNode->right = internalNodesQueue.front();
1016             prntNode->right->parent = prntNode;
1017             internalNodesQueue.pop();
1018         }
1019
1020         if ( lidx <= 0)
1021         {
1022             prntNode->left = cldNode = data->new_node( 0, 0, 0, 0 );
1023             *leafValsuesIt >> cldNode->value; leafValsuesIt--;
1024             cldNode->parent = prntNode;
1025         }
1026         else
1027         {
1028             prntNode->left = internalNodesQueue.front();
1029             prntNode->left->parent = prntNode;
1030             internalNodesQueue.pop();
1031         }
1032
1033         internalNodesQueue.push( prntNode );
1034     }
1035
1036     root = internalNodesQueue.front();
1037     internalNodesQueue.pop();
1038 }
1039
1040 void CvCascadeBoostTree::split_node_data( CvDTreeNode* node )
1041 {
1042     int n = node->sample_count, nl, nr, scount = data->sample_count;
1043     char* dir = (char*)data->direction->data.ptr;
1044     CvDTreeNode *left = 0, *right = 0;
1045     int* newIdx = data->split_buf->data.i;
1046     int newBufIdx = data->get_child_buf_idx( node );
1047     int workVarCount = data->get_work_var_count();
1048     CvMat* buf = data->buf;
1049     size_t length_buf_row = data->get_length_subbuf();
1050     cv::AutoBuffer<uchar> inn_buf(n*(3*sizeof(int)+sizeof(float)));
1051     int* tempBuf = (int*)(uchar*)inn_buf;
1052     bool splitInputData;
1053
1054     complete_node_dir(node);
1055
1056     for( int i = nl = nr = 0; i < n; i++ )
1057     {
1058         int d = dir[i];
1059         // initialize new indices for splitting ordered variables
1060         newIdx[i] = (nl & (d-1)) | (nr & -d); // d ? ri : li
1061         nr += d;
1062         nl += d^1;
1063     }
1064
1065     node->left = left = data->new_node( node, nl, newBufIdx, node->offset );
1066     node->right = right = data->new_node( node, nr, newBufIdx, node->offset + nl );
1067
1068     splitInputData = node->depth + 1 < data->params.max_depth &&
1069         (node->left->sample_count > data->params.min_sample_count ||
1070         node->right->sample_count > data->params.min_sample_count);
1071
1072     // split ordered variables, keep both halves sorted.
1073     for( int vi = 0; vi < ((CvCascadeBoostTrainData*)data)->numPrecalcIdx; vi++ )
1074     {
1075         int ci = data->get_var_type(vi);
1076         if( ci >= 0 || !splitInputData )
1077             continue;
1078
1079         int n1 = node->get_num_valid(vi);
1080         float *src_val_buf = (float*)(tempBuf + n);
1081         int *src_sorted_idx_buf = (int*)(src_val_buf + n);
1082         int *src_sample_idx_buf = src_sorted_idx_buf + n;
1083         const int* src_sorted_idx = 0;
1084         const float* src_val = 0;
1085         data->get_ord_var_data(node, vi, src_val_buf, src_sorted_idx_buf, &src_val, &src_sorted_idx, src_sample_idx_buf);
1086
1087         for(int i = 0; i < n; i++)
1088             tempBuf[i] = src_sorted_idx[i];
1089
1090         if (data->is_buf_16u)
1091         {
1092             ushort *ldst, *rdst;
1093             ldst = (ushort*)(buf->data.s + left->buf_idx*length_buf_row +
1094                 vi*scount + left->offset);
1095             rdst = (ushort*)(ldst + nl);
1096
1097             // split sorted
1098             for( int i = 0; i < n1; i++ )
1099             {
1100                 int idx = tempBuf[i];
1101                 int d = dir[idx];
1102                 idx = newIdx[idx];
1103                 if (d)
1104                 {
1105                     *rdst = (ushort)idx;
1106                     rdst++;
1107                 }
1108                 else
1109                 {
1110                     *ldst = (ushort)idx;
1111                     ldst++;
1112                 }
1113             }
1114             CV_Assert( n1 == n );
1115         }
1116         else
1117         {
1118             int *ldst, *rdst;
1119             ldst = buf->data.i + left->buf_idx*length_buf_row +
1120                 vi*scount + left->offset;
1121             rdst = buf->data.i + right->buf_idx*length_buf_row +
1122                 vi*scount + right->offset;
1123
1124             // split sorted
1125             for( int i = 0; i < n1; i++ )
1126             {
1127                 int idx = tempBuf[i];
1128                 int d = dir[idx];
1129                 idx = newIdx[idx];
1130                 if (d)
1131                 {
1132                     *rdst = idx;
1133                     rdst++;
1134                 }
1135                 else
1136                 {
1137                     *ldst = idx;
1138                     ldst++;
1139                 }
1140             }
1141             CV_Assert( n1 == n );
1142         }
1143     }
1144
1145     // split cv_labels using newIdx relocation table
1146     int *src_lbls_buf = tempBuf + n;
1147     const int* src_lbls = data->get_cv_labels(node, src_lbls_buf);
1148
1149     for(int i = 0; i < n; i++)
1150         tempBuf[i] = src_lbls[i];
1151
1152     if (data->is_buf_16u)
1153     {
1154         unsigned short *ldst = (unsigned short *)(buf->data.s + left->buf_idx*length_buf_row +
1155             (workVarCount-1)*scount + left->offset);
1156         unsigned short *rdst = (unsigned short *)(buf->data.s + right->buf_idx*length_buf_row +
1157             (workVarCount-1)*scount + right->offset);
1158
1159         for( int i = 0; i < n; i++ )
1160         {
1161             int idx = tempBuf[i];
1162             if (dir[i])
1163             {
1164                 *rdst = (unsigned short)idx;
1165                 rdst++;
1166             }
1167             else
1168             {
1169                 *ldst = (unsigned short)idx;
1170                 ldst++;
1171             }
1172         }
1173
1174     }
1175     else
1176     {
1177         int *ldst = buf->data.i + left->buf_idx*length_buf_row +
1178             (workVarCount-1)*scount + left->offset;
1179         int *rdst = buf->data.i + right->buf_idx*length_buf_row +
1180             (workVarCount-1)*scount + right->offset;
1181
1182         for( int i = 0; i < n; i++ )
1183         {
1184             int idx = tempBuf[i];
1185             if (dir[i])
1186             {
1187                 *rdst = idx;
1188                 rdst++;
1189             }
1190             else
1191             {
1192                 *ldst = idx;
1193                 ldst++;
1194             }
1195         }
1196     }
1197
1198     // split sample indices
1199     int *sampleIdx_src_buf = tempBuf + n;
1200     const int* sampleIdx_src = data->get_sample_indices(node, sampleIdx_src_buf);
1201
1202     for(int i = 0; i < n; i++)
1203         tempBuf[i] = sampleIdx_src[i];
1204
1205     if (data->is_buf_16u)
1206     {
1207         unsigned short* ldst = (unsigned short*)(buf->data.s + left->buf_idx*length_buf_row +
1208             workVarCount*scount + left->offset);
1209         unsigned short* rdst = (unsigned short*)(buf->data.s + right->buf_idx*length_buf_row +
1210             workVarCount*scount + right->offset);
1211         for (int i = 0; i < n; i++)
1212         {
1213             unsigned short idx = (unsigned short)tempBuf[i];
1214             if (dir[i])
1215             {
1216                 *rdst = idx;
1217                 rdst++;
1218             }
1219             else
1220             {
1221                 *ldst = idx;
1222                 ldst++;
1223             }
1224         }
1225     }
1226     else
1227     {
1228         int* ldst = buf->data.i + left->buf_idx*length_buf_row +
1229             workVarCount*scount + left->offset;
1230         int* rdst = buf->data.i + right->buf_idx*length_buf_row +
1231             workVarCount*scount + right->offset;
1232         for (int i = 0; i < n; i++)
1233         {
1234             int idx = tempBuf[i];
1235             if (dir[i])
1236             {
1237                 *rdst = idx;
1238                 rdst++;
1239             }
1240             else
1241             {
1242                 *ldst = idx;
1243                 ldst++;
1244             }
1245         }
1246     }
1247
1248     for( int vi = 0; vi < data->var_count; vi++ )
1249     {
1250         left->set_num_valid(vi, (int)(nl));
1251         right->set_num_valid(vi, (int)(nr));
1252     }
1253
1254     // deallocate the parent node data that is not needed anymore
1255     data->free_node_data(node);
1256 }
1257
1258 static void auxMarkFeaturesInMap( const CvDTreeNode* node, Mat& featureMap)
1259 {
1260     if ( node && node->split )
1261     {
1262         featureMap.ptr<int>(0)[node->split->var_idx] = 1;
1263         auxMarkFeaturesInMap( node->left, featureMap );
1264         auxMarkFeaturesInMap( node->right, featureMap );
1265     }
1266 }
1267
1268 void CvCascadeBoostTree::markFeaturesInMap( Mat& featureMap )
1269 {
1270     auxMarkFeaturesInMap( root, featureMap );
1271 }
1272
1273 //----------------------------------- CascadeBoost --------------------------------------
1274
1275 bool CvCascadeBoost::train( const CvFeatureEvaluator* _featureEvaluator,
1276                            int _numSamples,
1277                            int _precalcValBufSize, int _precalcIdxBufSize,
1278                            const CvCascadeBoostParams& _params )
1279 {
1280     bool isTrained = false;
1281     CV_Assert( !data );
1282     clear();
1283     data = new CvCascadeBoostTrainData( _featureEvaluator, _numSamples,
1284                                         _precalcValBufSize, _precalcIdxBufSize, _params );
1285     CvMemStorage *storage = cvCreateMemStorage();
1286     weak = cvCreateSeq( 0, sizeof(CvSeq), sizeof(CvBoostTree*), storage );
1287     storage = 0;
1288
1289     set_params( _params );
1290     if ( (_params.boost_type == LOGIT) || (_params.boost_type == GENTLE) )
1291         data->do_responses_copy();
1292
1293     update_weights( 0 );
1294
1295     cout << "+----+---------+---------+" << endl;
1296     cout << "|  N |    HR   |    FA   |" << endl;
1297     cout << "+----+---------+---------+" << endl;
1298
1299     do
1300     {
1301         CvCascadeBoostTree* tree = new CvCascadeBoostTree;
1302         if( !tree->train( data, subsample_mask, this ) )
1303         {
1304             delete tree;
1305             break;
1306         }
1307         cvSeqPush( weak, &tree );
1308         update_weights( tree );
1309         trim_weights();
1310         if( cvCountNonZero(subsample_mask) == 0 )
1311             break;
1312     }
1313     while( !isErrDesired() && (weak->total < params.weak_count) );
1314
1315     if(weak->total > 0)
1316     {
1317         data->is_classifier = true;
1318         data->free_train_data();
1319         isTrained = true;
1320     }
1321     else
1322         clear();
1323
1324     return isTrained;
1325 }
1326
1327 float CvCascadeBoost::predict( int sampleIdx, bool returnSum ) const
1328 {
1329     CV_Assert( weak );
1330     double sum = 0;
1331     CvSeqReader reader;
1332     cvStartReadSeq( weak, &reader );
1333     cvSetSeqReaderPos( &reader, 0 );
1334     for( int i = 0; i < weak->total; i++ )
1335     {
1336         CvBoostTree* wtree;
1337         CV_READ_SEQ_ELEM( wtree, reader );
1338         sum += ((CvCascadeBoostTree*)wtree)->predict(sampleIdx)->value;
1339     }
1340     if( !returnSum )
1341         sum = sum < threshold - CV_THRESHOLD_EPS ? 0.0 : 1.0;
1342     return (float)sum;
1343 }
1344
1345 bool CvCascadeBoost::set_params( const CvBoostParams& _params )
1346 {
1347     minHitRate = ((CvCascadeBoostParams&)_params).minHitRate;
1348     maxFalseAlarm = ((CvCascadeBoostParams&)_params).maxFalseAlarm;
1349     return ( ( minHitRate > 0 ) && ( minHitRate < 1) &&
1350         ( maxFalseAlarm > 0 ) && ( maxFalseAlarm < 1) &&
1351         CvBoost::set_params( _params ));
1352 }
1353
1354 void CvCascadeBoost::update_weights( CvBoostTree* tree )
1355 {
1356     int n = data->sample_count;
1357     double sumW = 0.;
1358     int step = 0;
1359     float* fdata = 0;
1360     int *sampleIdxBuf;
1361     const int* sampleIdx = 0;
1362     int inn_buf_size = ((params.boost_type == LOGIT) || (params.boost_type == GENTLE) ? n*sizeof(int) : 0) +
1363                        ( !tree ? n*sizeof(int) : 0 );
1364     cv::AutoBuffer<uchar> inn_buf(inn_buf_size);
1365     uchar* cur_inn_buf_pos = (uchar*)inn_buf;
1366     if ( (params.boost_type == LOGIT) || (params.boost_type == GENTLE) )
1367     {
1368         step = CV_IS_MAT_CONT(data->responses_copy->type) ?
1369             1 : data->responses_copy->step / CV_ELEM_SIZE(data->responses_copy->type);
1370         fdata = data->responses_copy->data.fl;
1371         sampleIdxBuf = (int*)cur_inn_buf_pos; cur_inn_buf_pos = (uchar*)(sampleIdxBuf + n);
1372         sampleIdx = data->get_sample_indices( data->data_root, sampleIdxBuf );
1373     }
1374     CvMat* buf = data->buf;
1375     size_t length_buf_row = data->get_length_subbuf();
1376     if( !tree ) // before training the first tree, initialize weights and other parameters
1377     {
1378         int* classLabelsBuf = (int*)cur_inn_buf_pos; cur_inn_buf_pos = (uchar*)(classLabelsBuf + n);
1379         const int* classLabels = data->get_class_labels(data->data_root, classLabelsBuf);
1380         // in case of logitboost and gentle adaboost each weak tree is a regression tree,
1381         // so we need to convert class labels to floating-point values
1382         double w0 = 1./n;
1383         double p[2] = { 1, 1 };
1384
1385         cvReleaseMat( &orig_response );
1386         cvReleaseMat( &sum_response );
1387         cvReleaseMat( &weak_eval );
1388         cvReleaseMat( &subsample_mask );
1389         cvReleaseMat( &weights );
1390
1391         orig_response = cvCreateMat( 1, n, CV_32S );
1392         weak_eval = cvCreateMat( 1, n, CV_64F );
1393         subsample_mask = cvCreateMat( 1, n, CV_8U );
1394         weights = cvCreateMat( 1, n, CV_64F );
1395         subtree_weights = cvCreateMat( 1, n + 2, CV_64F );
1396
1397         if (data->is_buf_16u)
1398         {
1399             unsigned short* labels = (unsigned short*)(buf->data.s + data->data_root->buf_idx*length_buf_row +
1400                 data->data_root->offset + (data->work_var_count-1)*data->sample_count);
1401             for( int i = 0; i < n; i++ )
1402             {
1403                 // save original categorical responses {0,1}, convert them to {-1,1}
1404                 orig_response->data.i[i] = classLabels[i]*2 - 1;
1405                 // make all the samples active at start.
1406                 // later, in trim_weights() deactivate/reactive again some, if need
1407                 subsample_mask->data.ptr[i] = (uchar)1;
1408                 // make all the initial weights the same.
1409                 weights->data.db[i] = w0*p[classLabels[i]];
1410                 // set the labels to find (from within weak tree learning proc)
1411                 // the particular sample weight, and where to store the response.
1412                 labels[i] = (unsigned short)i;
1413             }
1414         }
1415         else
1416         {
1417             int* labels = buf->data.i + data->data_root->buf_idx*length_buf_row +
1418                 data->data_root->offset + (data->work_var_count-1)*data->sample_count;
1419
1420             for( int i = 0; i < n; i++ )
1421             {
1422                 // save original categorical responses {0,1}, convert them to {-1,1}
1423                 orig_response->data.i[i] = classLabels[i]*2 - 1;
1424                 subsample_mask->data.ptr[i] = (uchar)1;
1425                 weights->data.db[i] = w0*p[classLabels[i]];
1426                 labels[i] = i;
1427             }
1428         }
1429
1430         if( params.boost_type == LOGIT )
1431         {
1432             sum_response = cvCreateMat( 1, n, CV_64F );
1433
1434             for( int i = 0; i < n; i++ )
1435             {
1436                 sum_response->data.db[i] = 0;
1437                 fdata[sampleIdx[i]*step] = orig_response->data.i[i] > 0 ? 2.f : -2.f;
1438             }
1439
1440             // in case of logitboost each weak tree is a regression tree.
1441             // the target function values are recalculated for each of the trees
1442             data->is_classifier = false;
1443         }
1444         else if( params.boost_type == GENTLE )
1445         {
1446             for( int i = 0; i < n; i++ )
1447                 fdata[sampleIdx[i]*step] = (float)orig_response->data.i[i];
1448
1449             data->is_classifier = false;
1450         }
1451     }
1452     else
1453     {
1454         // at this moment, for all the samples that participated in the training of the most
1455         // recent weak classifier we know the responses. For other samples we need to compute them
1456         if( have_subsample )
1457         {
1458             // invert the subsample mask
1459             cvXorS( subsample_mask, cvScalar(1.), subsample_mask );
1460
1461             // run tree through all the non-processed samples
1462             for( int i = 0; i < n; i++ )
1463                 if( subsample_mask->data.ptr[i] )
1464                 {
1465                     weak_eval->data.db[i] = ((CvCascadeBoostTree*)tree)->predict( i )->value;
1466                 }
1467         }
1468
1469         // now update weights and other parameters for each type of boosting
1470         if( params.boost_type == DISCRETE )
1471         {
1472             // Discrete AdaBoost:
1473             //   weak_eval[i] (=f(x_i)) is in {-1,1}
1474             //   err = sum(w_i*(f(x_i) != y_i))/sum(w_i)
1475             //   C = log((1-err)/err)
1476             //   w_i *= exp(C*(f(x_i) != y_i))
1477
1478             double C, err = 0.;
1479             double scale[] = { 1., 0. };
1480
1481             for( int i = 0; i < n; i++ )
1482             {
1483                 double w = weights->data.db[i];
1484                 sumW += w;
1485                 err += w*(weak_eval->data.db[i] != orig_response->data.i[i]);
1486             }
1487
1488             if( sumW != 0 )
1489                 err /= sumW;
1490             C = err = -logRatio( err );
1491             scale[1] = exp(err);
1492
1493             sumW = 0;
1494             for( int i = 0; i < n; i++ )
1495             {
1496                 double w = weights->data.db[i]*
1497                     scale[weak_eval->data.db[i] != orig_response->data.i[i]];
1498                 sumW += w;
1499                 weights->data.db[i] = w;
1500             }
1501
1502             tree->scale( C );
1503         }
1504         else if( params.boost_type == REAL )
1505         {
1506             // Real AdaBoost:
1507             //   weak_eval[i] = f(x_i) = 0.5*log(p(x_i)/(1-p(x_i))), p(x_i)=P(y=1|x_i)
1508             //   w_i *= exp(-y_i*f(x_i))
1509
1510             for( int i = 0; i < n; i++ )
1511                 weak_eval->data.db[i] *= -orig_response->data.i[i];
1512
1513             cvExp( weak_eval, weak_eval );
1514
1515             for( int i = 0; i < n; i++ )
1516             {
1517                 double w = weights->data.db[i]*weak_eval->data.db[i];
1518                 sumW += w;
1519                 weights->data.db[i] = w;
1520             }
1521         }
1522         else if( params.boost_type == LOGIT )
1523         {
1524             // LogitBoost:
1525             //   weak_eval[i] = f(x_i) in [-z_max,z_max]
1526             //   sum_response = F(x_i).
1527             //   F(x_i) += 0.5*f(x_i)
1528             //   p(x_i) = exp(F(x_i))/(exp(F(x_i)) + exp(-F(x_i))=1/(1+exp(-2*F(x_i)))
1529             //   reuse weak_eval: weak_eval[i] <- p(x_i)
1530             //   w_i = p(x_i)*1(1 - p(x_i))
1531             //   z_i = ((y_i+1)/2 - p(x_i))/(p(x_i)*(1 - p(x_i)))
1532             //   store z_i to the data->data_root as the new target responses
1533
1534             const double lbWeightThresh = FLT_EPSILON;
1535             const double lbZMax = 10.;
1536
1537             for( int i = 0; i < n; i++ )
1538             {
1539                 double s = sum_response->data.db[i] + 0.5*weak_eval->data.db[i];
1540                 sum_response->data.db[i] = s;
1541                 weak_eval->data.db[i] = -2*s;
1542             }
1543
1544             cvExp( weak_eval, weak_eval );
1545
1546             for( int i = 0; i < n; i++ )
1547             {
1548                 double p = 1./(1. + weak_eval->data.db[i]);
1549                 double w = p*(1 - p), z;
1550                 w = MAX( w, lbWeightThresh );
1551                 weights->data.db[i] = w;
1552                 sumW += w;
1553                 if( orig_response->data.i[i] > 0 )
1554                 {
1555                     z = 1./p;
1556                     fdata[sampleIdx[i]*step] = (float)min(z, lbZMax);
1557                 }
1558                 else
1559                 {
1560                     z = 1./(1-p);
1561                     fdata[sampleIdx[i]*step] = (float)-min(z, lbZMax);
1562                 }
1563             }
1564         }
1565         else
1566         {
1567             // Gentle AdaBoost:
1568             //   weak_eval[i] = f(x_i) in [-1,1]
1569             //   w_i *= exp(-y_i*f(x_i))
1570             assert( params.boost_type == GENTLE );
1571
1572             for( int i = 0; i < n; i++ )
1573                 weak_eval->data.db[i] *= -orig_response->data.i[i];
1574
1575             cvExp( weak_eval, weak_eval );
1576
1577             for( int i = 0; i < n; i++ )
1578             {
1579                 double w = weights->data.db[i] * weak_eval->data.db[i];
1580                 weights->data.db[i] = w;
1581                 sumW += w;
1582             }
1583         }
1584     }
1585
1586     // renormalize weights
1587     if( sumW > FLT_EPSILON )
1588     {
1589         sumW = 1./sumW;
1590         for( int i = 0; i < n; ++i )
1591             weights->data.db[i] *= sumW;
1592     }
1593 }
1594
1595 bool CvCascadeBoost::isErrDesired()
1596 {
1597     int sCount = data->sample_count,
1598         numPos = 0, numNeg = 0, numFalse = 0, numPosTrue = 0;
1599     vector<float> eval(sCount);
1600
1601     for( int i = 0; i < sCount; i++ )
1602         if( ((CvCascadeBoostTrainData*)data)->featureEvaluator->getCls( i ) == 1.0F )
1603             eval[numPos++] = predict( i, true );
1604     icvSortFlt( &eval[0], numPos, 0 );
1605     int thresholdIdx = (int)((1.0F - minHitRate) * numPos);
1606     threshold = eval[ thresholdIdx ];
1607     numPosTrue = numPos - thresholdIdx;
1608     for( int i = thresholdIdx - 1; i >= 0; i--)
1609         if ( abs( eval[i] - threshold) < FLT_EPSILON )
1610             numPosTrue++;
1611     float hitRate = ((float) numPosTrue) / ((float) numPos);
1612
1613     for( int i = 0; i < sCount; i++ )
1614     {
1615         if( ((CvCascadeBoostTrainData*)data)->featureEvaluator->getCls( i ) == 0.0F )
1616         {
1617             numNeg++;
1618             if( predict( i ) )
1619                 numFalse++;
1620         }
1621     }
1622     float falseAlarm = ((float) numFalse) / ((float) numNeg);
1623
1624     cout << "|"; cout.width(4); cout << right << weak->total;
1625     cout << "|"; cout.width(9); cout << right << hitRate;
1626     cout << "|"; cout.width(9); cout << right << falseAlarm;
1627     cout << "|" << endl;
1628     cout << "+----+---------+---------+" << endl;
1629
1630     return falseAlarm <= maxFalseAlarm;
1631 }
1632
1633 void CvCascadeBoost::write( FileStorage &fs, const Mat& featureMap ) const
1634 {
1635 //    char cmnt[30];
1636     CvCascadeBoostTree* weakTree;
1637     fs << CC_WEAK_COUNT << weak->total;
1638     fs << CC_STAGE_THRESHOLD << threshold;
1639     fs << CC_WEAK_CLASSIFIERS << "[";
1640     for( int wi = 0; wi < weak->total; wi++)
1641     {
1642         /*sprintf( cmnt, "tree %i", wi );
1643         cvWriteComment( fs, cmnt, 0 );*/
1644         weakTree = *((CvCascadeBoostTree**) cvGetSeqElem( weak, wi ));
1645         weakTree->write( fs, featureMap );
1646     }
1647     fs << "]";
1648 }
1649
1650 bool CvCascadeBoost::read( const FileNode &node,
1651                            const CvFeatureEvaluator* _featureEvaluator,
1652                            const CvCascadeBoostParams& _params )
1653 {
1654     CvMemStorage* storage;
1655     clear();
1656     data = new CvCascadeBoostTrainData( _featureEvaluator, _params );
1657     set_params( _params );
1658
1659     node[CC_STAGE_THRESHOLD] >> threshold;
1660     FileNode rnode = node[CC_WEAK_CLASSIFIERS];
1661
1662     storage = cvCreateMemStorage();
1663     weak = cvCreateSeq( 0, sizeof(CvSeq), sizeof(CvBoostTree*), storage );
1664     for( FileNodeIterator it = rnode.begin(); it != rnode.end(); it++ )
1665     {
1666         CvCascadeBoostTree* tree = new CvCascadeBoostTree();
1667         tree->read( *it, this, data );
1668         cvSeqPush( weak, &tree );
1669     }
1670     return true;
1671 }
1672
1673 void CvCascadeBoost::markUsedFeaturesInMap( Mat& featureMap )
1674 {
1675     for( int wi = 0; wi < weak->total; wi++ )
1676     {
1677         CvCascadeBoostTree* weakTree = *((CvCascadeBoostTree**) cvGetSeqElem( weak, wi ));
1678         weakTree->markFeaturesInMap( featureMap );
1679     }
1680 }