opencv/apps/traincascade/boost.cpp

   1 #include "boost.h"
   2 #include "cascadeclassifier.h"
   3 #include <queue>
   4
   5 #ifdef _OPENMP
   6 #include "omp.h"
   7 #endif
   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 //----------------------------- CascadeBoostParams -------------------------------------------------
  34
  35 CvCascadeBoostParams::CvCascadeBoostParams() : minHitRate( 0.995F), maxFalseAlarm( 0.5F )
  36 {
  37     boost_type = CvBoost::GENTLE;
  38     use_surrogates = use_1se_rule = truncate_pruned_tree = false;
  39 }
  40
  41 CvCascadeBoostParams::CvCascadeBoostParams( int _boostType,
  42         float _minHitRate, float _maxFalseAlarm,
  43         double _weightTrimRate, int _maxDepth, int _maxWeakCount ) :
  44     CvBoostParams( _boostType, _maxWeakCount, _weightTrimRate, _maxDepth, false, 0 )
  45 {
  46     boost_type = CvBoost::GENTLE;
  47     minHitRate = _minHitRate;
  48     maxFalseAlarm = _maxFalseAlarm;
  49     use_surrogates = use_1se_rule = truncate_pruned_tree = false;
  50 }
  51
  52 void CvCascadeBoostParams::write( FileStorage &fs ) const
  53 {
  54     String boostTypeStr = boost_type == CvBoost::DISCRETE ? CC_DISCRETE_BOOST :
  55                           boost_type == CvBoost::REAL ? CC_REAL_BOOST :
  56                           boost_type == CvBoost::LOGIT ? CC_LOGIT_BOOST :
  57                           boost_type == CvBoost::GENTLE ? CC_GENTLE_BOOST : String();
  58     CV_Assert( !boostTypeStr.empty() );
  59     fs << CC_BOOST_TYPE << boostTypeStr;
  60     fs << CC_MINHITRATE << minHitRate;
  61     fs << CC_MAXFALSEALARM << maxFalseAlarm;
  62     fs << CC_TRIM_RATE << weight_trim_rate;
  63     fs << CC_MAX_DEPTH << max_depth;
  64     fs << CC_WEAK_COUNT << weak_count;
  65 }
  66
  67 bool CvCascadeBoostParams::read( const FileNode &node )
  68 {
  69     String boostTypeStr;
  70     FileNode rnode = node[CC_BOOST_TYPE];
  71     rnode >> boostTypeStr;
  72     boost_type = !boostTypeStr.compare( CC_DISCRETE_BOOST ) ? CvBoost::DISCRETE :
  73                  !boostTypeStr.compare( CC_REAL_BOOST ) ? CvBoost::REAL :
  74                  !boostTypeStr.compare( CC_LOGIT_BOOST ) ? CvBoost::LOGIT :
  75                  !boostTypeStr.compare( CC_GENTLE_BOOST ) ? CvBoost::GENTLE : -1;
  76     if (boost_type == -1)
  77         CV_Error( CV_StsBadArg, "unsupported Boost type" );
  78     node[CC_MINHITRATE] >> minHitRate;
  79     node[CC_MAXFALSEALARM] >> maxFalseAlarm;
  80     node[CC_TRIM_RATE] >> weight_trim_rate ;
  81     node[CC_MAX_DEPTH] >> max_depth ;
  82     node[CC_WEAK_COUNT] >> weak_count ;
  83     if ( minHitRate <= 0 || minHitRate > 1 ||
  84          maxFalseAlarm <= 0 || maxFalseAlarm > 1 ||
  85          weight_trim_rate <= 0 || weight_trim_rate > 1 ||
  86          max_depth <= 0 || weak_count <= 0 )
  87         CV_Error( CV_StsBadArg, "bad parameters range");
  88     return true;
  89 }
  90
  91 void CvCascadeBoostParams::printDefaults() const
  92 {
  93     cout << "--boostParams--" << endl;
  94     cout << "  [-bt <{" << CC_DISCRETE_BOOST << ", "
  95                         << CC_REAL_BOOST << ", "
  96                         << CC_LOGIT_BOOST ", "
  97                         << CC_GENTLE_BOOST << "(default)}>]" << endl;
  98     cout << "  [-minHitRate <min_hit_rate> = " << minHitRate << ">]" << endl;
  99     cout << "  [-maxFalseAlarmRate <max_false_alarm_rate = " << maxFalseAlarm << ">]" << endl;
 100     cout << "  [-weightTrimRate <weight_trim_rate = " << weight_trim_rate << ">]" << endl;
 101     cout << "  [-maxDepth <max_depth_of_weak_tree = " << max_depth << ">]" << endl;
 102     cout << "  [-maxWeakCount <max_weak_tree_count = " << weak_count << ">]" << endl;
 103 }
 104
 105 void CvCascadeBoostParams::printAttrs() const
 106 {
 107     String boostTypeStr = boost_type == CvBoost::DISCRETE ? CC_DISCRETE_BOOST :
 108                           boost_type == CvBoost::REAL ? CC_REAL_BOOST :
 109                           boost_type == CvBoost::LOGIT  ? CC_LOGIT_BOOST :
 110                           boost_type == CvBoost::GENTLE ? CC_GENTLE_BOOST : String();
 111     CV_Assert( !boostTypeStr.empty() );
 112     cout << "boostType: " << boostTypeStr << endl;
 113     cout << "minHitRate: " << minHitRate << endl;
 114     cout << "maxFalseAlarmRate: " <<  maxFalseAlarm << endl;
 115     cout << "weightTrimRate: " << weight_trim_rate << endl;
 116     cout << "maxTreeDepth: " << max_depth << endl;
 117     cout << "maxWeakCount: " << weak_count << endl;
 118 }
 119
 120 bool CvCascadeBoostParams::scanAttr( const String prmName, const String val)
 121 {
 122     bool res = true;
 123
 124     if( !prmName.compare( "-bt" ) )
 125     {
 126         boost_type = !val.compare( CC_DISCRETE_BOOST ) ? CvBoost::DISCRETE :
 127                      !val.compare( CC_REAL_BOOST ) ? CvBoost::REAL :
 128                      !val.compare( CC_LOGIT_BOOST ) ? CvBoost::LOGIT :
 129                      !val.compare( CC_GENTLE_BOOST ) ? CvBoost::GENTLE : -1;
 130         if (boost_type == -1)
 131             res = false;
 132     }
 133     else if( !prmName.compare( "-minHitRate" ) )
 134     {
 135         minHitRate = (float) atof( val.c_str() );
 136     }
 137     else if( !prmName.compare( "-maxFalseAlarmRate" ) )
 138     {
 139         weight_trim_rate = (float) atof( val.c_str() );
 140     }
 141     else if( !prmName.compare( "-weightTrimRate" ) )
 142     {
 143         weight_trim_rate = (float) atof( val.c_str() );
 144     }
 145     else if( !prmName.compare( "-maxDepth" ) )
 146     {
 147         max_depth = atoi( val.c_str() );
 148     }
 149     else if( !prmName.compare( "-maxWeakCount" ) )
 150     {
 151         weak_count = atoi( val.c_str() );
 152     }
 153     else
 154         res = false;
 155
 156     return res;
 157 }
 158
 159 //---------------------------- CascadeBoostTrainData -----------------------------
 160
 161 CvCascadeBoostTrainData::CvCascadeBoostTrainData( const CvFeatureEvaluator* _featureEvaluator,
 162                                                   const CvDTreeParams& _params )
 163 {
 164     is_classifier = true;
 165     var_all = var_count = (int)_featureEvaluator->getNumFeatures();
 166
 167     featureEvaluator = _featureEvaluator;
 168     shared = true;
 169     set_params( _params );
 170     max_c_count = MAX( 2, featureEvaluator->getMaxCatCount() );
 171     var_type = cvCreateMat( 1, var_count + 2, CV_32SC1 );
 172     if ( featureEvaluator->getMaxCatCount() > 0 )
 173     {
 174         numPrecalcIdx = 0;
 175         cat_var_count = var_count;
 176         ord_var_count = 0;
 177         for( int vi = 0; vi < var_count; vi++ )
 178         {
 179             var_type->data.i[vi] = vi;
 180         }
 181     }
 182     else
 183     {
 184         cat_var_count = 0;
 185         ord_var_count = var_count;
 186         for( int vi = 1; vi <= var_count; vi++ )
 187         {
 188             var_type->data.i[vi-1] = -vi;
 189         }
 190     }
 191     var_type->data.i[var_count] = cat_var_count;
 192     var_type->data.i[var_count+1] = cat_var_count+1;
 193
 194     int maxSplitSize = cvAlign(sizeof(CvDTreeSplit) + (MAX(0,max_c_count - 33)/32)*sizeof(int),sizeof(void*));
 195     int treeBlockSize = MAX((int)sizeof(CvDTreeNode)*8, maxSplitSize);
 196     treeBlockSize = MAX(treeBlockSize + BlockSizeDelta, MinBlockSize);
 197     tree_storage = cvCreateMemStorage( treeBlockSize );
 198     node_heap = cvCreateSet( 0, sizeof(node_heap[0]), sizeof(CvDTreeNode), tree_storage );
 199     split_heap = cvCreateSet( 0, sizeof(split_heap[0]), maxSplitSize, tree_storage );
 200 }
 201
 202 CvCascadeBoostTrainData::CvCascadeBoostTrainData( const CvFeatureEvaluator* _featureEvaluator,
 203                                                  int _numSamples,
 204                                                  int _precalcValBufSize, int _precalcIdxBufSize,
 205                                                  const CvDTreeParams& _params )
 206 {
 207     setData( _featureEvaluator, _numSamples, _precalcValBufSize, _precalcIdxBufSize, _params );
 208 }
 209
 210 void CvCascadeBoostTrainData::setData( const CvFeatureEvaluator* _featureEvaluator,
 211                                       int _numSamples,
 212                                       int _precalcValBufSize, int _precalcIdxBufSize,
 213                                                                           const CvDTreeParams& _params )
 214 {
 215     int* idst = 0;
 216     unsigned short* udst = 0;
 217
 218     clear();
 219     shared = true;
 220     have_labels = true;
 221     have_priors = false;
 222     is_classifier = true;
 223
 224     rng = cvRNG(-1);
 225
 226     set_params( _params );
 227
 228     CV_Assert( _featureEvaluator );
 229     featureEvaluator = _featureEvaluator;
 230
 231     max_c_count = MAX( 2, featureEvaluator->getMaxCatCount() );
 232     _resp = featureEvaluator->getCls();
 233     responses = &_resp;
 234     // TODO: check responses: elements must be 0 or 1
 235
 236         if( _precalcValBufSize < 0 || _precalcIdxBufSize < 0)
 237         CV_Error( CV_StsOutOfRange, "_numPrecalcVal and _numPrecalcIdx must be positive or 0" );
 238
 239         var_count = var_all = featureEvaluator->getNumFeatures();
 240     sample_count = _numSamples;
 241
 242     is_buf_16u = false;
 243     if (sample_count < 65536)
 244         is_buf_16u = true;
 245
 246         numPrecalcVal = min( (_precalcValBufSize*1048576) / int(sizeof(float)*sample_count), var_count );
 247     numPrecalcIdx = min( (_precalcIdxBufSize*1048576) /
 248                 int((is_buf_16u ? sizeof(unsigned short) : sizeof (int))*sample_count), var_count );
 249
 250     int maxNumThreads = 1;
 251 #ifdef _OPENMP
 252     maxNumThreads = omp_get_num_procs();
 253 #endif
 254     valCache.create( max(numPrecalcVal, maxNumThreads), sample_count, CV_32FC1 );
 255     var_type = cvCreateMat( 1, var_count + 2, CV_32SC1 );
 256
 257     if ( featureEvaluator->getMaxCatCount() > 0 )
 258     {
 259         numPrecalcIdx = 0;
 260         cat_var_count = var_count;
 261         ord_var_count = 0;
 262         for( int vi = 0; vi < var_count; vi++ )
 263         {
 264             var_type->data.i[vi] = vi;
 265         }
 266     }
 267     else
 268     {
 269         cat_var_count = 0;
 270         ord_var_count = var_count;
 271         for( int vi = 1; vi <= var_count; vi++ )
 272         {
 273             var_type->data.i[vi-1] = -vi;
 274         }
 275     }
 276     var_type->data.i[var_count] = cat_var_count;
 277     var_type->data.i[var_count+1] = cat_var_count+1;
 278     work_var_count = ( cat_var_count ? 0 : numPrecalcIdx ) + 1;
 279     buf_size = (work_var_count + 1) * sample_count;
 280     buf_count = 2;
 281
 282     if ( is_buf_16u )
 283         buf = cvCreateMat( buf_count, buf_size, CV_16UC1 );
 284     else
 285         buf = cvCreateMat( buf_count, buf_size, CV_32SC1 );
 286
 287     cat_count = cvCreateMat( 1, cat_var_count + 1, CV_32SC1 );
 288     pred_float_buf.resize(maxNumThreads);
 289     pred_int_buf.resize(maxNumThreads);
 290     resp_float_buf.resize(maxNumThreads);
 291     resp_int_buf.resize(maxNumThreads);
 292     cv_lables_buf.resize(maxNumThreads);
 293     sample_idx_buf.resize(maxNumThreads);
 294     for( int ti = 0; ti < maxNumThreads; ti++ )
 295     {
 296         pred_float_buf[ti].resize(sample_count);
 297         pred_int_buf[ti].resize(sample_count);
 298         resp_float_buf[ti].resize(sample_count);
 299         resp_int_buf[ti].resize(sample_count);
 300         cv_lables_buf[ti].resize(sample_count);
 301         sample_idx_buf[ti].resize(sample_count);
 302     }
 303
 304     // precalculate valCache and set indices in buf
 305         precalculate();
 306
 307     // now calculate the maximum size of split,
 308     // create memory storage that will keep nodes and splits of the decision tree
 309     // allocate root node and the buffer for the whole training data
 310     int maxSplitSize = cvAlign(sizeof(CvDTreeSplit) +
 311         (MAX(0,sample_count - 33)/32)*sizeof(int),sizeof(void*));
 312     int treeBlockSize = MAX((int)sizeof(CvDTreeNode)*8, maxSplitSize);
 313     treeBlockSize = MAX(treeBlockSize + BlockSizeDelta, MinBlockSize);
 314     tree_storage = cvCreateMemStorage( treeBlockSize );
 315     node_heap = cvCreateSet( 0, sizeof(*node_heap), sizeof(CvDTreeNode), tree_storage );
 316
 317     int nvSize = var_count*sizeof(int);
 318     nvSize = cvAlign(MAX( nvSize, (int)sizeof(CvSetElem) ), sizeof(void*));
 319     int tempBlockSize = nvSize;
 320     tempBlockSize = MAX( tempBlockSize + BlockSizeDelta, MinBlockSize );
 321     temp_storage = cvCreateMemStorage( tempBlockSize );
 322     nv_heap = cvCreateSet( 0, sizeof(*nv_heap), nvSize, temp_storage );
 323
 324     data_root = new_node( 0, sample_count, 0, 0 );
 325
 326     // set sample labels
 327     if (is_buf_16u)
 328         udst = (unsigned short*)(buf->data.s + work_var_count*sample_count);
 329     else
 330         idst = buf->data.i + work_var_count*sample_count;
 331
 332     for (int si = 0; si < sample_count; si++)
 333     {
 334         if (udst)
 335             udst[si] = (unsigned short)si;
 336         else
 337             idst[si] = si;
 338     }
 339     for( int vi = 0; vi < var_count; vi++ )
 340         data_root->set_num_valid(vi, sample_count);
 341     for( int vi = 0; vi < cat_var_count; vi++ )
 342         cat_count->data.i[vi] = max_c_count;
 343
 344     cat_count->data.i[cat_var_count] = 2;
 345
 346     maxSplitSize = cvAlign(sizeof(CvDTreeSplit) +
 347         (MAX(0,max_c_count - 33)/32)*sizeof(int),sizeof(void*));
 348     split_heap = cvCreateSet( 0, sizeof(*split_heap), maxSplitSize, tree_storage );
 349
 350     priors = cvCreateMat( 1, get_num_classes(), CV_64F );
 351     cvSet(priors, cvScalar(1));
 352     priors_mult = cvCloneMat( priors );
 353     counts = cvCreateMat( 1, get_num_classes(), CV_32SC1 );
 354     direction = cvCreateMat( 1, sample_count, CV_8UC1 );
 355     split_buf = cvCreateMat( 1, sample_count, CV_32SC1 );
 356 }
 357
 358 void CvCascadeBoostTrainData::free_train_data()
 359 {
 360     CvDTreeTrainData::free_train_data();
 361     valCache.release();
 362 }
 363
 364 void CvCascadeBoostTrainData::get_class_labels( CvDTreeNode* n, int* labelsBuf, const int** labels )
 365 {
 366     int nodeSampleCount = n->sample_count;
 367     int* sampleIndicesBuf = get_sample_idx_buf();
 368     const int* sampleIndices = 0;
 369     int rStep = CV_IS_MAT_CONT( responses->type ) ? 1 : responses->step / CV_ELEM_SIZE( responses->type );
 370
 371     get_sample_indices(n, sampleIndicesBuf, &sampleIndices);
 372
 373     for( int si = 0; si < nodeSampleCount; si++ )
 374     {
 375         int sidx = sampleIndices[si];
 376         labelsBuf[si] = (int)responses->data.fl[sidx*rStep];
 377     }
 378     *labels = labelsBuf;
 379 }
 380
 381 void CvCascadeBoostTrainData::get_sample_indices( CvDTreeNode* n, int* indicesBuf, const int** indices )
 382 {
 383     CvDTreeTrainData::get_cat_var_data( n, get_work_var_count(), indicesBuf, indices );
 384 }
 385
 386 void CvCascadeBoostTrainData::get_cv_labels( CvDTreeNode* n, int* labels_buf, const int** labels )
 387 {
 388     CvDTreeTrainData::get_cat_var_data( n, get_work_var_count()- 1, labels_buf, labels );
 389 }
 390
 391 int CvCascadeBoostTrainData::get_ord_var_data( CvDTreeNode* n, int vi, float* ordValuesBuf, int* indicesBuf,
 392         const float** ordValues, const int** indices )
 393 {
 394     int nodeSampleCount = n->sample_count;
 395     int* sampleIndicesBuf = get_sample_idx_buf();
 396     const int* sampleIndices = 0;
 397     get_sample_indices(n, sampleIndicesBuf, &sampleIndices);
 398
 399         if ( vi < numPrecalcIdx )
 400         {
 401                 if( !is_buf_16u )
 402                 *indices = buf->data.i + n->buf_idx*buf->cols + vi*sample_count + n->offset;
 403             else
 404         {
 405                 const unsigned short* shortIndices = (const unsigned short*)(buf->data.s + n->buf_idx*buf->cols +
 406                                                       vi*sample_count + n->offset );
 407                 for( int i = 0; i < nodeSampleCount; i++ )
 408                     indicesBuf[i] = shortIndices[i];
 409                 *indices = indicesBuf;
 410             }
 411
 412                 if ( vi < numPrecalcVal )
 413                 {
 414                         for( int i = 0; i < nodeSampleCount; i++ )
 415                 {
 416                     int idx = (*indices)[i];
 417                     idx = sampleIndices[idx];
 418                     ordValuesBuf[i] =  valCache.at<float>( vi, idx);
 419                 }
 420                 }
 421                 else
 422                 {
 423                         for( int i = 0; i < nodeSampleCount; i++ )
 424                 {
 425                     int idx = (*indices)[i];
 426                     idx = sampleIndices[idx];
 427                                 ordValuesBuf[i] = (*featureEvaluator)( vi, idx);
 428                         }
 429                 }
 430     }
 431         else // vi >= numPrecalcIdx
 432         {
 433                 // use sample_indices as temporary buffer for values
 434                 if ( vi < numPrecalcVal )
 435                 {
 436                         for( int i = 0; i < nodeSampleCount; i++ )
 437                 {
 438                                 indicesBuf[i] = i;
 439                     ((float*)sampleIndices)[i] = valCache.at<float>( vi, sampleIndices[i] );
 440                 }
 441                 }
 442                 else
 443                 {
 444             for( int i = 0; i < nodeSampleCount; i++ )
 445                 {
 446                                 indicesBuf[i] = i;
 447                                 ((float*)sampleIndices)[i] = (*featureEvaluator)( vi, sampleIndices[i]);
 448                         }
 449                 }
 450                 icvSortIntAux( indicesBuf, sample_count, (float *)sampleIndices );
 451                 for( int i = 0; i < nodeSampleCount; i++ )
 452                 ordValuesBuf[i] = ((float*)sampleIndices)[indicesBuf[i]];
 453                 *indices = indicesBuf;
 454         }
 455
 456     *ordValues = ordValuesBuf;
 457     return 0;
 458 }
 459
 460 int CvCascadeBoostTrainData::get_cat_var_data( CvDTreeNode* n, int vi, int* catValuesBuf, const int** catValues )
 461 {
 462     int nodeSampleCount = n->sample_count;
 463     int* sampleIndicesBuf = get_sample_idx_buf();
 464     const int* sampleIndices = 0;
 465     get_sample_indices(n, sampleIndicesBuf, &sampleIndices);
 466
 467     if ( vi < numPrecalcVal )
 468         {
 469                 for( int i = 0; i < nodeSampleCount; i++ )
 470             catValuesBuf[i] = (int) valCache.at<float>( vi, sampleIndices[i]);
 471         }
 472         else
 473         {
 474                 for( int i = 0; i < nodeSampleCount; i++ )
 475                         catValuesBuf[i] = (int)(*featureEvaluator)( vi, sampleIndices[i] );
 476         }
 477     *catValues = catValuesBuf;
 478
 479     return 0;
 480 }
 481
 482 float CvCascadeBoostTrainData::getVarValue( int vi, int si )
 483 {
 484     if ( vi < numPrecalcVal && !valCache.empty() )
 485             return valCache.at<float>( vi, si );
 486         return (*featureEvaluator)( vi, si );
 487 }
 488
 489 void CvCascadeBoostTrainData::precalculate()
 490 {
 491     int minNum = MIN( numPrecalcVal, numPrecalcIdx);
 492     unsigned short* udst = (unsigned short*)buf->data.s;
 493     int* idst = buf->data.i;
 494
 495     CV_DbgAssert( !valCache.empty() );
 496     double proctime = -TIME( 0 );
 497
 498 #ifdef _OPENMP
 499     int maxNumThreads = omp_get_num_procs();
 500 #pragma omp parallel for num_threads(maxNumThreads) schedule(dynamic)
 501 #endif
 502     for ( int fi = numPrecalcVal; fi < numPrecalcIdx; fi++)
 503     {
 504         int threadIdx = 0;
 505 #ifdef _OPENMP
 506                 threadIdx = omp_get_thread_num();
 507 #endif
 508         for( int si = 0; si < sample_count; si++ )
 509         {
 510             valCache.ptr<float>(threadIdx)[si] = (*featureEvaluator)( fi, si );
 511             if ( is_buf_16u )
 512                 *(udst + fi*sample_count + si) = (unsigned short)si;
 513             else
 514                 *(idst + fi*sample_count + si) = si;
 515         }
 516         if ( is_buf_16u )
 517             icvSortUShAux( udst + fi*sample_count, sample_count, valCache.ptr<float>(threadIdx) );
 518         else
 519             icvSortIntAux( idst + fi*sample_count, sample_count, valCache.ptr<float>(threadIdx) );
 520     }
 521 #ifdef _OPENMP
 522 #pragma omp parallel for num_threads(maxNumThreads) schedule(dynamic)
 523 #endif
 524     for ( int fi = 0; fi < minNum; fi++)
 525     {
 526         for( int si = 0; si < sample_count; si++ )
 527         {
 528             valCache.ptr<float>(fi)[si] = (*featureEvaluator)( fi, si );
 529             if ( is_buf_16u )
 530                 *(udst + fi*sample_count + si) = (unsigned short)si;
 531             else
 532                 *(idst + fi*sample_count + si) = si;
 533         }
 534         if ( is_buf_16u )
 535             icvSortUShAux( udst + fi*sample_count, sample_count, valCache.ptr<float>(fi) );
 536         else
 537             icvSortIntAux( idst + fi*sample_count, sample_count, valCache.ptr<float>(fi) );
 538     }
 539 #ifdef _OPENMP
 540 #pragma omp parallel for num_threads(maxNumThreads) schedule(dynamic)
 541 #endif
 542     for ( int fi = minNum; fi < numPrecalcVal; fi++)
 543         for( int si = 0; si < sample_count; si++ )
 544             valCache.ptr<float>(fi)[si] = (*featureEvaluator)( fi, si );
 545
 546     cout << "Precalculation time: " << (proctime + TIME( 0 )) << endl;
 547 }
 548
 549 //-------------------------------- CascadeBoostTree ----------------------------------------
 550
 551 CvDTreeNode* CvCascadeBoostTree::predict( int sampleIdx ) const
 552 {
 553     CvDTreeNode* node = root;
 554     if( !node )
 555         CV_Error( CV_StsError, "The tree has not been trained yet" );
 556
 557     if ( ((CvCascadeBoostTrainData*)data)->featureEvaluator->getMaxCatCount() == 0 ) // ordered
 558     {
 559         while( node->left )
 560         {
 561             CvDTreeSplit* split = node->split;
 562             float val = ((CvCascadeBoostTrainData*)data)->getVarValue( split->var_idx, sampleIdx );
 563             node = val <= split->ord.c ? node->left : node->right;
 564         }
 565     }
 566     else // categorical
 567     {
 568         while( node->left )
 569         {
 570             CvDTreeSplit* split = node->split;
 571             int c = (int)((CvCascadeBoostTrainData*)data)->getVarValue( split->var_idx, sampleIdx );
 572             node = CV_DTREE_CAT_DIR(c, split->subset) < 0 ? node->left : node->right;
 573         }
 574     }
 575     return node;
 576 }
 577
 578 void CvCascadeBoostTree::write( FileStorage &fs, const Mat& featureMap )
 579 {
 580     int maxCatCount = ((CvCascadeBoostTrainData*)data)->featureEvaluator->getMaxCatCount();
 581     int subsetN = (maxCatCount + 31)/32;
 582     queue<CvDTreeNode*> internalNodesQueue;
 583     int size = (int)pow( 2.f, (float)ensemble->get_params().max_depth);
 584     Ptr<float> leafVals = new float[size];
 585     int leafValIdx = 0;
 586     int internalNodeIdx = 1;
 587     CvDTreeNode* tempNode;
 588
 589     CV_DbgAssert( root );
 590     internalNodesQueue.push( root );
 591
 592     fs << "{";
 593     fs << CC_INTERNAL_NODES << "[:";
 594     while (!internalNodesQueue.empty())
 595     {
 596         tempNode = internalNodesQueue.front();
 597         CV_Assert( tempNode->left );
 598         if ( !tempNode->left->left && !tempNode->left->right) // left node is leaf
 599         {
 600             leafVals[-leafValIdx] = (float)tempNode->left->value;
 601             fs << leafValIdx-- ;
 602         }
 603         else
 604         {
 605             internalNodesQueue.push( tempNode->left );
 606             fs << internalNodeIdx++;
 607         }
 608         CV_Assert( tempNode->right );
 609         if ( !tempNode->right->left && !tempNode->right->right) // right node is leaf
 610         {
 611             leafVals[-leafValIdx] = (float)tempNode->right->value;
 612             fs << leafValIdx--;
 613         }
 614         else
 615         {
 616             internalNodesQueue.push( tempNode->right );
 617             fs << internalNodeIdx++;
 618         }
 619         int fidx = tempNode->split->var_idx;
 620         fidx = featureMap.empty() ? fidx : featureMap.at<int>(0, fidx);
 621         fs << fidx;
 622         if ( !maxCatCount )
 623             fs << tempNode->split->ord.c;
 624         else
 625             for( int i = 0; i < subsetN; i++ )
 626                 fs << tempNode->split->subset[i];
 627         internalNodesQueue.pop();
 628     }
 629     fs << "]"; // CC_INTERNAL_NODES
 630
 631     fs << CC_LEAF_VALUES << "[:";
 632     for (int ni = 0; ni < -leafValIdx; ni++)
 633         fs << leafVals[ni];
 634     fs << "]"; // CC_LEAF_VALUES
 635     fs << "}";
 636 }
 637
 638 void CvCascadeBoostTree::read( const FileNode &node, CvBoost* _ensemble,
 639                                 CvDTreeTrainData* _data )
 640 {
 641     int maxCatCount = ((CvCascadeBoostTrainData*)_data)->featureEvaluator->getMaxCatCount();
 642     int subsetN = (maxCatCount + 31)/32;
 643     int step = 3 + ( maxCatCount>0 ? subsetN : 1 );
 644
 645     queue<CvDTreeNode*> internalNodesQueue;
 646     FileNodeIterator internalNodesIt, leafValsuesIt;
 647     CvDTreeNode* prntNode, *cldNode;
 648
 649     clear();
 650     data = _data;
 651     ensemble = _ensemble;
 652     pruned_tree_idx = 0;
 653
 654     // read tree nodes
 655     FileNode rnode = node[CC_INTERNAL_NODES];
 656     internalNodesIt = rnode.end();
 657     leafValsuesIt = node[CC_LEAF_VALUES].end();
 658     internalNodesIt--; leafValsuesIt--;
 659     for( size_t i = 0; i < rnode.size()/step; i++ )
 660     {
 661         prntNode = data->new_node( 0, 0, 0, 0 );
 662         if ( maxCatCount > 0 )
 663         {
 664             prntNode->split = data->new_split_cat( 0, 0 );
 665             for( int j = subsetN-1; j>=0; j--)
 666             {
 667                 *internalNodesIt >> prntNode->split->subset[j]; internalNodesIt--;
 668             }
 669         }
 670         else
 671         {
 672             float split_value;
 673             *internalNodesIt >> split_value; internalNodesIt--;
 674             prntNode->split = data->new_split_ord( 0, split_value, 0, 0, 0);
 675         }
 676         *internalNodesIt >> prntNode->split->var_idx; internalNodesIt--;
 677         int ridx, lidx;
 678         *internalNodesIt >> ridx; internalNodesIt--;
 679         *internalNodesIt >> lidx;internalNodesIt--;
 680         if ( ridx <= 0)
 681         {
 682             prntNode->right = cldNode = data->new_node( 0, 0, 0, 0 );
 683             *leafValsuesIt >> cldNode->value; leafValsuesIt--;
 684             cldNode->parent = prntNode;
 685         }
 686         else
 687         {
 688             prntNode->right = internalNodesQueue.front();
 689             prntNode->right->parent = prntNode;
 690             internalNodesQueue.pop();
 691         }
 692
 693         if ( lidx <= 0)
 694         {
 695             prntNode->left = cldNode = data->new_node( 0, 0, 0, 0 );
 696             *leafValsuesIt >> cldNode->value; leafValsuesIt--;
 697             cldNode->parent = prntNode;
 698         }
 699         else
 700         {
 701             prntNode->left = internalNodesQueue.front();
 702             prntNode->left->parent = prntNode;
 703             internalNodesQueue.pop();
 704         }
 705
 706         internalNodesQueue.push( prntNode );
 707     }
 708
 709     root = internalNodesQueue.front();
 710     internalNodesQueue.pop();
 711 }
 712
 713 void CvCascadeBoostTree::split_node_data( CvDTreeNode* node )
 714 {
 715     int n = node->sample_count, nl, nr, scount = data->sample_count;
 716     char* dir = (char*)data->direction->data.ptr;
 717     CvDTreeNode *left = 0, *right = 0;
 718     int* newIdx = data->split_buf->data.i;
 719     int newBufIdx = data->get_child_buf_idx( node );
 720     int workVarCount = data->get_work_var_count();
 721     CvMat* buf = data->buf;
 722     int* tempBuf = (int*)cvStackAlloc(n*sizeof(tempBuf[0]));
 723     bool splitInputData;
 724
 725     complete_node_dir(node);
 726
 727     for( int i = nl = nr = 0; i < n; i++ )
 728     {
 729         int d = dir[i];
 730         // initialize new indices for splitting ordered variables
 731         newIdx[i] = (nl & (d-1)) | (nr & -d); // d ? ri : li
 732         nr += d;
 733         nl += d^1;
 734     }
 735
 736     node->left = left = data->new_node( node, nl, newBufIdx, node->offset );
 737     node->right = right = data->new_node( node, nr, newBufIdx, node->offset + nl );
 738
 739     splitInputData = node->depth + 1 < data->params.max_depth &&
 740         (node->left->sample_count > data->params.min_sample_count ||
 741         node->right->sample_count > data->params.min_sample_count);
 742
 743         // split ordered variables, keep both halves sorted.
 744     for( int vi = 0; vi < ((CvCascadeBoostTrainData*)data)->numPrecalcIdx; vi++ )
 745     {
 746         int ci = data->get_var_type(vi);
 747         int n1 = node->get_num_valid(vi);
 748         int *src_idx_buf = data->get_pred_int_buf();
 749         const int* src_idx = 0;
 750         float *src_val_buf = data->get_pred_float_buf();
 751         const float* src_val = 0;
 752
 753         if( ci >= 0 || !splitInputData )
 754             continue;
 755
 756         data->get_ord_var_data(node, vi, src_val_buf, src_idx_buf, &src_val, &src_idx);
 757
 758         for(int i = 0; i < n; i++)
 759             tempBuf[i] = src_idx[i];
 760
 761         if (data->is_buf_16u)
 762         {
 763             unsigned short *ldst, *rdst, *ldst0, *rdst0;
 764             ldst0 = ldst = (unsigned short*)(buf->data.s + left->buf_idx*buf->cols +
 765                 vi*scount + left->offset);
 766             rdst0 = rdst = (unsigned short*)(ldst + nl);
 767
 768             // split sorted
 769             for( int i = 0; i < n1; i++ )
 770             {
 771                 int idx = tempBuf[i];
 772                 int d = dir[idx];
 773                 idx = newIdx[idx];
 774                 if (d)
 775                 {
 776                     *rdst = (unsigned short)idx;
 777                     rdst++;
 778                 }
 779                 else
 780                 {
 781                     *ldst = (unsigned short)idx;
 782                     ldst++;
 783                 }
 784             }
 785             assert( n1 == n);
 786
 787             left->set_num_valid(vi, (int)(ldst - ldst0));
 788             right->set_num_valid(vi, (int)(rdst - rdst0));
 789         }
 790         else
 791         {
 792             int *ldst0, *ldst, *rdst0, *rdst;
 793             ldst0 = ldst = buf->data.i + left->buf_idx*buf->cols +
 794                 vi*scount + left->offset;
 795             rdst0 = rdst = buf->data.i + right->buf_idx*buf->cols +
 796                 vi*scount + right->offset;
 797
 798             // split sorted
 799             for( int i = 0; i < n1; i++ )
 800             {
 801                 int idx = tempBuf[i];
 802                 int d = dir[idx];
 803                 idx = newIdx[idx];
 804                 if (d)
 805                 {
 806                     *rdst = idx;
 807                     rdst++;
 808                 }
 809                 else
 810                 {
 811                     *ldst = idx;
 812                     ldst++;
 813                 }
 814             }
 815
 816             left->set_num_valid(vi, (int)(ldst - ldst0));
 817             right->set_num_valid(vi, (int)(rdst - rdst0));
 818             CV_Assert( n1 == n);
 819         }
 820     }
 821
 822     // split cv_labels using newIdx relocation table
 823     int *src_lbls_buf = data->get_pred_int_buf();
 824     const int* src_lbls = 0;
 825     data->get_cv_labels(node, src_lbls_buf, &src_lbls);
 826
 827     for(int i = 0; i < n; i++)
 828         tempBuf[i] = src_lbls[i];
 829
 830     if (data->is_buf_16u)
 831     {
 832         unsigned short *ldst = (unsigned short *)(buf->data.s + left->buf_idx*buf->cols +
 833             (workVarCount-1)*scount + left->offset);
 834         unsigned short *rdst = (unsigned short *)(buf->data.s + right->buf_idx*buf->cols +
 835             (workVarCount-1)*scount + right->offset);
 836
 837         for( int i = 0; i < n; i++ )
 838         {
 839             int idx = tempBuf[i];
 840             if (dir[i])
 841             {
 842                 *rdst = (unsigned short)idx;
 843                 rdst++;
 844             }
 845             else
 846             {
 847                 *ldst = (unsigned short)idx;
 848                 ldst++;
 849             }
 850         }
 851
 852     }
 853     else
 854     {
 855         int *ldst = buf->data.i + left->buf_idx*buf->cols +
 856             (workVarCount-1)*scount + left->offset;
 857         int *rdst = buf->data.i + right->buf_idx*buf->cols +
 858             (workVarCount-1)*scount + right->offset;
 859
 860         for( int i = 0; i < n; i++ )
 861         {
 862             int idx = tempBuf[i];
 863             if (dir[i])
 864             {
 865                 *rdst = idx;
 866                 rdst++;
 867             }
 868             else
 869             {
 870                 *ldst = idx;
 871                 ldst++;
 872             }
 873         }
 874     }
 875     for( int vi = 0; vi < data->var_count; vi++ )
 876     {
 877         left->set_num_valid(vi, (int)(nl));
 878         right->set_num_valid(vi, (int)(nr));
 879     }
 880
 881     // split sample indices
 882     int *sampleIdx_src_buf = data->get_sample_idx_buf();
 883     const int* sampleIdx_src = 0;
 884     data->get_sample_indices(node, sampleIdx_src_buf, &sampleIdx_src);
 885
 886     for(int i = 0; i < n; i++)
 887         tempBuf[i] = sampleIdx_src[i];
 888
 889     if (data->is_buf_16u)
 890     {
 891         unsigned short* ldst = (unsigned short*)(buf->data.s + left->buf_idx*buf->cols +
 892             workVarCount*scount + left->offset);
 893         unsigned short* rdst = (unsigned short*)(buf->data.s + right->buf_idx*buf->cols +
 894             workVarCount*scount + right->offset);
 895         for (int i = 0; i < n; i++)
 896         {
 897             unsigned short idx = (unsigned short)tempBuf[i];
 898             if (dir[i])
 899             {
 900                 *rdst = idx;
 901                 rdst++;
 902             }
 903             else
 904             {
 905                 *ldst = idx;
 906                 ldst++;
 907             }
 908         }
 909     }
 910     else
 911     {
 912         int* ldst = buf->data.i + left->buf_idx*buf->cols +
 913             workVarCount*scount + left->offset;
 914         int* rdst = buf->data.i + right->buf_idx*buf->cols +
 915             workVarCount*scount + right->offset;
 916         for (int i = 0; i < n; i++)
 917         {
 918             int idx = tempBuf[i];
 919             if (dir[i])
 920             {
 921                 *rdst = idx;
 922                 rdst++;
 923             }
 924             else
 925             {
 926                 *ldst = idx;
 927                 ldst++;
 928             }
 929         }
 930     }
 931
 932     // deallocate the parent node data that is not needed anymore
 933     data->free_node_data(node);
 934 }
 935
 936 void auxMarkFeaturesInMap( const CvDTreeNode* node, Mat& featureMap)
 937 {
 938     if ( node && node->split )
 939     {
 940         featureMap.ptr<int>(0)[node->split->var_idx] = 1;
 941         auxMarkFeaturesInMap( node->left, featureMap );
 942         auxMarkFeaturesInMap( node->right, featureMap );
 943     }
 944 }
 945
 946 void CvCascadeBoostTree::markFeaturesInMap( Mat& featureMap )
 947 {
 948     auxMarkFeaturesInMap( root, featureMap );
 949 }
 950
 951 //----------------------------------- CascadeBoost --------------------------------------
 952
 953 bool CvCascadeBoost::train( const CvFeatureEvaluator* _featureEvaluator,
 954                            int _numSamples,
 955                            int _precalcValBufSize, int _precalcIdxBufSize,
 956                            const CvCascadeBoostParams& _params )
 957 {
 958     CV_Assert( !data );
 959     clear();
 960     data = new CvCascadeBoostTrainData( _featureEvaluator, _numSamples,
 961                                         _precalcValBufSize, _precalcIdxBufSize, _params );
 962     CvMemStorage *storage = cvCreateMemStorage();
 963     weak = cvCreateSeq( 0, sizeof(CvSeq), sizeof(CvBoostTree*), storage );
 964     storage = 0;
 965
 966     set_params( _params );
 967     if ( (_params.boost_type == LOGIT) || (_params.boost_type == GENTLE) )
 968         data->do_responses_copy();
 969
 970     update_weights( 0 );
 971
 972     cout << "+----+---------+---------+" << endl;
 973     cout << "|  N |    HR   |    FA   |" << endl;
 974     cout << "+----+---------+---------+" << endl;
 975
 976     do
 977     {
 978         CvCascadeBoostTree* tree = new CvCascadeBoostTree;
 979         if( !tree->train( data, subsample_mask, this ) )
 980         {
 981             // TODO: may be should finish the loop (!!!)
 982             assert(0);
 983             delete tree;
 984             continue;
 985         }
 986         cvSeqPush( weak, &tree );
 987         update_weights( tree );
 988         trim_weights();
 989     }
 990     while( !isErrDesired() && (weak->total < params.weak_count) );
 991
 992     data->is_classifier = true;
 993     data->free_train_data();
 994     return true;
 995 }
 996
 997 float CvCascadeBoost::predict( int sampleIdx, bool returnSum ) const
 998 {
 999     CV_Assert( weak );
1000     double sum = 0;
1001     CvSeqReader reader;
1002     cvStartReadSeq( weak, &reader );
1003     cvSetSeqReaderPos( &reader, 0 );
1004     for( int i = 0; i < weak->total; i++ )
1005     {
1006         CvBoostTree* wtree;
1007         CV_READ_SEQ_ELEM( wtree, reader );
1008         sum += ((CvCascadeBoostTree*)wtree)->predict(sampleIdx)->value;
1009     }
1010     if( !returnSum )
1011         sum = sum < threshold - CV_THRESHOLD_EPS ? 0.0 : 1.0;
1012     return (float)sum;
1013 }
1014
1015 bool CvCascadeBoost::set_params( const CvBoostParams& _params )
1016 {
1017     minHitRate = ((CvCascadeBoostParams&)_params).minHitRate;
1018     maxFalseAlarm = ((CvCascadeBoostParams&)_params).maxFalseAlarm;
1019     return ( ( minHitRate > 0 ) && ( minHitRate < 1) &&
1020         ( maxFalseAlarm > 0 ) && ( maxFalseAlarm < 1) &&
1021         CvBoost::set_params( _params ));
1022 }
1023
1024 void CvCascadeBoost::update_weights( CvBoostTree* tree )
1025 {
1026     int n = data->sample_count;
1027     double sumW = 0.;
1028     int step = 0;
1029     float* fdata = 0;
1030     int *sampleIdxBuf;
1031     const int* sampleIdx = 0;
1032     if ( (params.boost_type == LOGIT) || (params.boost_type == GENTLE) )
1033     {
1034         step = CV_IS_MAT_CONT(data->responses_copy->type) ?
1035             1 : data->responses_copy->step / CV_ELEM_SIZE(data->responses_copy->type);
1036         fdata = data->responses_copy->data.fl;
1037         sampleIdxBuf = (int*)cvStackAlloc(data->sample_count*sizeof(sampleIdxBuf[0]));
1038         data->get_sample_indices( data->data_root, sampleIdxBuf, &sampleIdx );
1039     }
1040     CvMat* buf = data->buf;
1041     if( !tree ) // before training the first tree, initialize weights and other parameters
1042     {
1043         int* classLabelsBuf = data->get_resp_int_buf();
1044         const int* classLabels = 0;
1045         data->get_class_labels(data->data_root, classLabelsBuf, &classLabels);
1046         // in case of logitboost and gentle adaboost each weak tree is a regression tree,
1047         // so we need to convert class labels to floating-point values
1048         float* responses_buf = data->get_resp_float_buf();
1049         const float* responses = 0;
1050         data->get_ord_responses(data->data_root, responses_buf, &responses);
1051
1052         double w0 = 1./n;
1053         double p[2] = { 1, 1 };
1054
1055         cvReleaseMat( &orig_response );
1056         cvReleaseMat( &sum_response );
1057         cvReleaseMat( &weak_eval );
1058         cvReleaseMat( &subsample_mask );
1059         cvReleaseMat( &weights );
1060
1061         orig_response = cvCreateMat( 1, n, CV_32S );
1062         weak_eval = cvCreateMat( 1, n, CV_64F );
1063         subsample_mask = cvCreateMat( 1, n, CV_8U );
1064         weights = cvCreateMat( 1, n, CV_64F );
1065         subtree_weights = cvCreateMat( 1, n + 2, CV_64F );
1066
1067         if (data->is_buf_16u)
1068         {
1069             unsigned short* labels = (unsigned short*)(buf->data.s + data->data_root->buf_idx*buf->cols +
1070                 data->data_root->offset + (data->work_var_count-1)*data->sample_count);
1071             for( int i = 0; i < n; i++ )
1072             {
1073                 // save original categorical responses {0,1}, convert them to {-1,1}
1074                 orig_response->data.i[i] = classLabels[i]*2 - 1;
1075                 // make all the samples active at start.
1076                 // later, in trim_weights() deactivate/reactive again some, if need
1077                 subsample_mask->data.ptr[i] = (uchar)1;
1078                 // make all the initial weights the same.
1079                 weights->data.db[i] = w0*p[classLabels[i]];
1080                 // set the labels to find (from within weak tree learning proc)
1081                 // the particular sample weight, and where to store the response.
1082                 labels[i] = (unsigned short)i;
1083             }
1084         }
1085         else
1086         {
1087             int* labels = buf->data.i + data->data_root->buf_idx*buf->cols +
1088                 data->data_root->offset + (data->work_var_count-1)*data->sample_count;
1089
1090             for( int i = 0; i < n; i++ )
1091             {
1092                 // save original categorical responses {0,1}, convert them to {-1,1}
1093                 orig_response->data.i[i] = classLabels[i]*2 - 1;
1094                 subsample_mask->data.ptr[i] = (uchar)1;
1095                 weights->data.db[i] = w0*p[classLabels[i]];
1096                 labels[i] = i;
1097             }
1098         }
1099
1100         if( params.boost_type == LOGIT )
1101         {
1102             sum_response = cvCreateMat( 1, n, CV_64F );
1103
1104             for( int i = 0; i < n; i++ )
1105             {
1106                 sum_response->data.db[i] = 0;
1107                 fdata[sampleIdx[i]*step] = orig_response->data.i[i] > 0 ? 2.f : -2.f;
1108             }
1109
1110             // in case of logitboost each weak tree is a regression tree.
1111             // the target function values are recalculated for each of the trees
1112             data->is_classifier = false;
1113         }
1114         else if( params.boost_type == GENTLE )
1115         {
1116             for( int i = 0; i < n; i++ )
1117                 fdata[sampleIdx[i]*step] = (float)orig_response->data.i[i];
1118
1119             data->is_classifier = false;
1120         }
1121     }
1122     else
1123     {
1124         // at this moment, for all the samples that participated in the training of the most
1125         // recent weak classifier we know the responses. For other samples we need to compute them
1126         if( have_subsample )
1127         {
1128             // invert the subsample mask
1129             cvXorS( subsample_mask, cvScalar(1.), subsample_mask );
1130
1131             // run tree through all the non-processed samples
1132             for( int i = 0; i < n; i++ )
1133                 if( subsample_mask->data.ptr[i] )
1134                 {
1135                     weak_eval->data.db[i] = ((CvCascadeBoostTree*)tree)->predict( i )->value;
1136                 }
1137         }
1138
1139         // now update weights and other parameters for each type of boosting
1140         if( params.boost_type == DISCRETE )
1141         {
1142             // Discrete AdaBoost:
1143             //   weak_eval[i] (=f(x_i)) is in {-1,1}
1144             //   err = sum(w_i*(f(x_i) != y_i))/sum(w_i)
1145             //   C = log((1-err)/err)
1146             //   w_i *= exp(C*(f(x_i) != y_i))
1147
1148             double C, err = 0.;
1149             double scale[] = { 1., 0. };
1150
1151             for( int i = 0; i < n; i++ )
1152             {
1153                 double w = weights->data.db[i];
1154                 sumW += w;
1155                 err += w*(weak_eval->data.db[i] != orig_response->data.i[i]);
1156             }
1157
1158             if( sumW != 0 )
1159                 err /= sumW;
1160             C = err = -logRatio( err );
1161             scale[1] = exp(err);
1162
1163             sumW = 0;
1164             for( int i = 0; i < n; i++ )
1165             {
1166                 double w = weights->data.db[i]*
1167                     scale[weak_eval->data.db[i] != orig_response->data.i[i]];
1168                 sumW += w;
1169                 weights->data.db[i] = w;
1170             }
1171
1172             tree->scale( C );
1173         }
1174         else if( params.boost_type == REAL )
1175         {
1176             // Real AdaBoost:
1177             //   weak_eval[i] = f(x_i) = 0.5*log(p(x_i)/(1-p(x_i))), p(x_i)=P(y=1|x_i)
1178             //   w_i *= exp(-y_i*f(x_i))
1179
1180             for( int i = 0; i < n; i++ )
1181                 weak_eval->data.db[i] *= -orig_response->data.i[i];
1182
1183             cvExp( weak_eval, weak_eval );
1184
1185             for( int i = 0; i < n; i++ )
1186             {
1187                 double w = weights->data.db[i]*weak_eval->data.db[i];
1188                 sumW += w;
1189                 weights->data.db[i] = w;
1190             }
1191         }
1192         else if( params.boost_type == LOGIT )
1193         {
1194             // LogitBoost:
1195             //   weak_eval[i] = f(x_i) in [-z_max,z_max]
1196             //   sum_response = F(x_i).
1197             //   F(x_i) += 0.5*f(x_i)
1198             //   p(x_i) = exp(F(x_i))/(exp(F(x_i)) + exp(-F(x_i))=1/(1+exp(-2*F(x_i)))
1199             //   reuse weak_eval: weak_eval[i] <- p(x_i)
1200             //   w_i = p(x_i)*1(1 - p(x_i))
1201             //   z_i = ((y_i+1)/2 - p(x_i))/(p(x_i)*(1 - p(x_i)))
1202             //   store z_i to the data->data_root as the new target responses
1203
1204             const double lbWeightThresh = FLT_EPSILON;
1205             const double lbZMax = 10.;
1206             float* responsesBuf = data->get_resp_float_buf();
1207             const float* responses = 0;
1208             data->get_ord_responses(data->data_root, responsesBuf, &responses);
1209
1210             for( int i = 0; i < n; i++ )
1211             {
1212                 double s = sum_response->data.db[i] + 0.5*weak_eval->data.db[i];
1213                 sum_response->data.db[i] = s;
1214                 weak_eval->data.db[i] = -2*s;
1215             }
1216
1217             cvExp( weak_eval, weak_eval );
1218
1219             for( int i = 0; i < n; i++ )
1220             {
1221                 double p = 1./(1. + weak_eval->data.db[i]);
1222                 double w = p*(1 - p), z;
1223                 w = MAX( w, lbWeightThresh );
1224                 weights->data.db[i] = w;
1225                 sumW += w;
1226                 if( orig_response->data.i[i] > 0 )
1227                 {
1228                     z = 1./p;
1229                     fdata[sampleIdx[i]*step] = (float)min(z, lbZMax);
1230                 }
1231                 else
1232                 {
1233                     z = 1./(1-p);
1234                     fdata[sampleIdx[i]*step] = (float)-min(z, lbZMax);
1235                 }
1236             }
1237         }
1238         else
1239         {
1240             // Gentle AdaBoost:
1241             //   weak_eval[i] = f(x_i) in [-1,1]
1242             //   w_i *= exp(-y_i*f(x_i))
1243             assert( params.boost_type == GENTLE );
1244
1245             for( int i = 0; i < n; i++ )
1246                 weak_eval->data.db[i] *= -orig_response->data.i[i];
1247
1248             cvExp( weak_eval, weak_eval );
1249
1250             for( int i = 0; i < n; i++ )
1251             {
1252                 double w = weights->data.db[i] * weak_eval->data.db[i];
1253                 weights->data.db[i] = w;
1254                 sumW += w;
1255             }
1256         }
1257     }
1258
1259     // renormalize weights
1260     if( sumW > FLT_EPSILON )
1261     {
1262         sumW = 1./sumW;
1263         for( int i = 0; i < n; ++i )
1264             weights->data.db[i] *= sumW;
1265     }
1266 }
1267
1268 bool CvCascadeBoost::isErrDesired()
1269 {
1270     int sCount = data->sample_count,
1271         numPos = 0, numNeg = 0, numFalse = 0, numPosTrue = 0;
1272     float* eval = (float*) cvStackAlloc( sizeof(eval[0]) * sCount );
1273
1274     for( int i = 0; i < sCount; i++ )
1275         if( ((CvCascadeBoostTrainData*)data)->featureEvaluator->getCls( i ) == 1.0F )
1276             eval[numPos++] = predict( i, true );
1277     icvSortFlt( eval, numPos, 0 );
1278     int thresholdIdx = (int)((1.0F - minHitRate) * numPos);
1279     threshold = eval[ thresholdIdx ];
1280     numPosTrue = numPos - thresholdIdx;
1281     for( int i = thresholdIdx - 1; i >= 0; i--)
1282         if ( abs( eval[i] - threshold) < FLT_EPSILON )
1283             numPosTrue++;
1284     float hitRate = ((float) numPosTrue) / ((float) numPos);
1285
1286     for( int i = 0; i < sCount; i++ )
1287     {
1288         if( ((CvCascadeBoostTrainData*)data)->featureEvaluator->getCls( i ) == 0.0F )
1289         {
1290             numNeg++;
1291             if( predict( i ) )
1292                 numFalse++;
1293         }
1294     }
1295     float falseAlarm = ((float) numFalse) / ((float) numNeg);
1296
1297     cout << "|"; cout.width(4); cout << right << weak->total;
1298     cout << "|"; cout.width(9); cout << right << hitRate;
1299     cout << "|"; cout.width(9); cout << right << falseAlarm;
1300     cout << "|" << endl;
1301     cout << "+----+---------+---------+" << endl;
1302
1303     return falseAlarm <= maxFalseAlarm;
1304 }
1305
1306 void CvCascadeBoost::write( FileStorage &fs, const Mat& featureMap ) const
1307 {
1308 //    char cmnt[30];
1309     CvCascadeBoostTree* weakTree;
1310     fs << CC_WEAK_COUNT << weak->total;
1311     fs << CC_STAGE_THRESHOLD << threshold;
1312     fs << CC_WEAK_CLASSIFIERS << "[";
1313     for( int wi = 0; wi < weak->total; wi++)
1314     {
1315         /*sprintf( cmnt, "tree %i", wi );
1316         cvWriteComment( fs, cmnt, 0 );*/
1317         weakTree = *((CvCascadeBoostTree**) cvGetSeqElem( weak, wi ));
1318         weakTree->write( fs, featureMap );
1319     }
1320     fs << "]";
1321 }
1322
1323 bool CvCascadeBoost::read( const FileNode &node,
1324                            const CvFeatureEvaluator* _featureEvaluator,
1325                            const CvCascadeBoostParams& _params )
1326 {
1327     CvMemStorage* storage;
1328     clear();
1329     data = new CvCascadeBoostTrainData( _featureEvaluator, _params );
1330     set_params( _params );
1331
1332     node[CC_STAGE_THRESHOLD] >> threshold;
1333     FileNode rnode = node[CC_WEAK_CLASSIFIERS];
1334
1335     storage = cvCreateMemStorage();
1336     weak = cvCreateSeq( 0, sizeof(CvSeq), sizeof(CvBoostTree*), storage );
1337     for( FileNodeIterator it = rnode.begin(); it != rnode.end(); it++ )
1338     {
1339         CvCascadeBoostTree* tree = new CvCascadeBoostTree();
1340         tree->read( *it, this, data );
1341         cvSeqPush( weak, &tree );
1342     }
1343     return true;
1344 }
1345
1346 void CvCascadeBoost::markUsedFeaturesInMap( Mat& featureMap )
1347 {
1348     for( int wi = 0; wi < weak->total; wi++ )
1349     {
1350         CvCascadeBoostTree* weakTree = *((CvCascadeBoostTree**) cvGetSeqElem( weak, wi ));
1351         weakTree->markFeaturesInMap( featureMap );
1352     }
1353 }