OpenCV 臉部偵測自學(6)opencv_traincascade如何訓練強弱分類器

在（3）中把opencv_traincascade在使用LBP特徵的時候的訓練準備工作的代碼總結了下。下面開始硬著頭皮看訓練裡面的部分了。介於這部分實在是沒怎麼找到網上介紹的文章（為啥呢？？？）所以總結的大部分內容是自己猜測的。以後再回頭慢慢完善。

接著上次結束的部分，訓練一個強分類器的代碼是：

bool CvCascadeBoost::train( const CvFeatureEvaluator* _featureEvaluator,//包含了sum，tilted，特徵的位置等資訊                           int _numSamples,                           int _precalcValBufSize, int _precalcIdxBufSize,                           const CvCascadeBoostParams& _params ){    CV_Assert( !data );    clear();    data = new CvCascadeBoostTrainData( _featureEvaluator, _numSamples,                                        _precalcValBufSize, _precalcIdxBufSize, _params );//目前暫且知道這裡是調用preCalculate計算所有的特徵值    CvMemStorage *storage = cvCreateMemStorage();    weak = cvCreateSeq( 0, sizeof(CvSeq), sizeof(CvBoostTree*), storage );    storage = 0;    set_params( _params );    if ( (_params.boost_type == LOGIT) || (_params.boost_type == GENTLE) )        data->do_responses_copy();    update_weights( 0 );    cout << "+----+---------+---------+" << endl;    cout << "|  N |    HR   |    FA   |" << endl;    cout << "+----+---------+---------+" << endl;    do    {        CvCascadeBoostTree* tree = new CvCascadeBoostTree;        if( !tree->train( data, subsample_mask, this ) )//應該是訓練一個弱分類器tree        {            delete tree;            break;        }        cvSeqPush( weak, &tree );//把弱分類器添加到強分類器裡面        update_weights( tree );//AdaBoost裡面更新weight部分。        trim_weights();        if( cvCountNonZero(subsample_mask) == 0 )            break;    }    while( !isErrDesired() && (weak->total < params.weak_count) );    data->is_classifier = true;    data->free_train_data();    return true;}

這裡面主要分為

  調用preCalculate計算所有樣本的所有的特徵值

      訓練一個弱分類器

      根據AdaBoost裡面的步驟把樣本的weight都更新一遍

  
主要看如何訓練一個弱分類器。因為在(4)中總結output的xml檔案結構的時候發現每一個node裡面都包含了8個很大的數，後來在predict函數裡看到這8個數是屬於subset，在使用xml檔案的時候用裡面每一個feature算出來的特徵值轉換為二進位後是256個bit，然後跟這8個大數的最後32個位進行比較，如果有至少一個bit都是1的話那麼就輸出一個值，反之另一個值。在這需要指出的是，如果是haar-like feature的話用的不是subset而是threshold，後面比較後就一樣了。而且haar-like基本上是一個feature一個弱分類器，而這裡subset裡可以有多個1出現，所以我理解著就是指有多個feature組成的subset。

而這些subset是如何求出來的呢，我跟蹤程式總結如下。（裡面用到好多opencv裡面ml的資料結構，到現在感覺還是沒有入門啊）

boost.cpp裡：

boolCvBoostTree::train( CvDTreeTrainData* _train_data,                    const CvMat* _subsample_idx, CvBoost* _ensemble ){    clear();    ensemble = _ensemble;    data = _train_data;    data->shared = true;    return do_train( _subsample_idx );//去CvDTree裡訓練弱分類器去}

然後在tree.cpp裡

bool CvDTree::do_train( const CvMat* _subsample_idx ){    bool result = false;    CV_FUNCNAME( "CvDTree::do_train" );    __BEGIN__;    root = data->subsample_data( _subsample_idx );    CV_CALL( try_split_node(root));//在這把root的value啊，爹媽孩子都算出來的，其中跟分類器的子集（threshold）有關的split->subset就在此計算的。    if( root->split )    {        CV_Assert( root->left );        CV_Assert( root->right );        if( data->params.cv_folds > 0 )            CV_CALL( prune_cv() );        if( !data->shared )            data->free_train_data();        result = true;    }    __END__;    return result;}

在這裡面主要包括：

    找root（回頭再研究）

   分裂root。 try_split_node(root)，就是在這把我想知道的那些數給算出來的，在opencv文檔上看到這麼一段話應該是對此的描述“決策樹是從根結點遞迴構造的。用所有的訓練資料（特徵向量和對應的響應）來在根結點處進行分裂。在每個結點處，最佳化準則（比如最優分裂）是基於一些基本原則來確定的（比如ML中的“純度purity”原則被用來進行分類，方差之和用來進行迴歸）。”

  prune_cv() 估計是剪枝（回頭再研究）

然後再研究分裂root部分。
tree.cpp中

void CvDTree::try_split_node( CvDTreeNode* node ){    CvDTreeSplit* best_split = 0;    int i, n = node->sample_count, vi;    bool can_split = true;    double quality_scale;    calc_node_value( node );//計算node->value    if( node->sample_count <= data->params.min_sample_count ||        node->depth >= data->params.max_depth )        can_split = false;    if( can_split && data->is_classifier )    {        // check if we have a "pure" node,        // we assume that cls_count is filled by calc_node_value()        int* cls_count = data->counts->data.i;        int nz = 0, m = data->get_num_classes();        for( i = 0; i < m; i++ )            nz += cls_count[i] != 0;        if( nz == 1 ) // there is only one class            can_split = false;    }    else if( can_split )    {        if( sqrt(node->node_risk)/n < data->params.regression_accuracy )            can_split = false;    }    if( can_split )    {        best_split = find_best_split(node);//可能是在這計算split->subset[]        // TODO: check the split quality ...        node->split = best_split;    }    if( !can_split || !best_split )    {        data->free_node_data(node);        return;    }    quality_scale = calc_node_dir( node );    if( data->params.use_surrogates )    {        // find all the surrogate splits        // and sort them by their similarity to the primary one        for( vi = 0; vi < data->var_count; vi++ )        {            CvDTreeSplit* split;            int ci = data->get_var_type(vi);            if( vi == best_split->var_idx )                continue;            if( ci >= 0 )                split = find_surrogate_split_cat( node, vi );            else                split = find_surrogate_split_ord( node, vi );            if( split )            {                // insert the split                CvDTreeSplit* prev_split = node->split;                split->quality = (float)(split->quality*quality_scale);                while( prev_split->next &&                       prev_split->next->quality > split->quality )                    prev_split = prev_split->next;                split->next = prev_split->next;                prev_split->next = split;            }        }    }    split_node_data( node );    try_split_node( node->left );///////    try_split_node( node->right );///////*演算法迴歸的繼續分裂左右孩子結點。在以下情況下演算法可能會在某個結點停止（i.e. stop splitting the node further）：樹的深度達到了指定的最大值在該結點訓練樣本的數目少於指定值，比如，沒有統計意義上的集合來進一步進行結點分裂了。在該結點所有的樣本屬於同一類（或者，如果是迴歸的話，變化已經非常小了）跟隨機播放相比，能選擇到的最好的分裂已經基本沒有什麼有意義的改進了。*/}

我要找的計算在best_split = find_best_split(node)中。

CvDTreeSplit* CvDTree::find_best_split( CvDTreeNode* node ){    DTreeBestSplitFinder finder( this, node );    cv::parallel_reduce(cv::BlockedRange(0, data->var_count), finder);//調用本cpp檔案中的DTreeBestSplitFinder::operator()    CvDTreeSplit *bestSplit = 0;    if( finder.bestSplit->quality > 0 )    {        bestSplit = data->new_split_cat( 0, -1.0f );        memcpy( bestSplit, finder.bestSplit, finder.splitSize );    }    return bestSplit;}

void DTreeBestSplitFinder::operator()(const BlockedRange& range){    int vi, vi1 = range.begin(), vi2 = range.end();    int n = node->sample_count;    CvDTreeTrainData* data = tree->get_data();    AutoBuffer<uchar> inn_buf(2*n*(sizeof(int) + sizeof(float)));    for( vi = vi1; vi < vi2; vi++ )//{0,8464}，就是一共有多少個feature，這個屬於外面的迴圈，對於每一個feature    {        CvDTreeSplit *res;        int ci = data->get_var_type(vi);        if( node->get_num_valid(vi) <= 1 )            continue;        if( data->is_classifier )        {            if( ci >= 0 )                res = tree->find_split_cat_class( node, vi, bestSplit->quality, split, (uchar*)inn_buf );            else                res = tree->find_split_ord_class( node, vi, bestSplit->quality, split, (uchar*)inn_buf );        }        else        {            if( ci >= 0 )res = tree->find_split_cat_reg( node, vi, bestSplit->quality, split, (uchar*)inn_buf );//計算split 調用的是CvDTreeSplit* CvBoostTree::find_split_cat_reg            else                res = tree->find_split_ord_reg( node, vi, bestSplit->quality, split, (uchar*)inn_buf );        }        if( res && bestSplit->quality < split->quality )                memcpy( (CvDTreeSplit*)bestSplit, (CvDTreeSplit*)split, splitSize );//更新bestSplit選取那個最好的feature    }}

LBP是屬於category的，所以調用res = tree->find_split_cat_reg( node, vi, bestSplit->quality, split, (uchar*)inn_buf )

boost.cpp中：

CvDTreeSplit*CvBoostTree::find_split_cat_reg( CvDTreeNode* node, int vi, float init_quality, CvDTreeSplit* _split, uchar* _ext_buf )//cat好像是categary的意思，reg是regression的意思。{    const double* weights = ensemble->get_subtree_weights()->data.db;    int ci = data->get_var_type(vi);    int n = node->sample_count;//正負樣本個數    int mi = data->cat_count->data.i[ci];    int base_size = (2*mi+3)*sizeof(double) + mi*sizeof(double*);    cv::AutoBuffer<uchar> inn_buf(base_size);    if( !_ext_buf )        inn_buf.allocate(base_size + n*(2*sizeof(int) + sizeof(float)));    uchar* base_buf = (uchar*)inn_buf;    uchar* ext_buf = _ext_buf ? _ext_buf : base_buf + base_size;    int* cat_labels_buf = (int*)ext_buf;    const int* cat_labels = data->get_cat_var_data(node, vi, cat_labels_buf);//返回的是1 x n的矩陣，裡面是所有樣本的第vi個feature的特徵值    float* responses_buf = (float*)(cat_labels_buf + n);    int* sample_indices_buf = (int*)(responses_buf + n);    const float* responses = data->get_ord_responses(node, responses_buf, sample_indices_buf);//response就是類別{1,-1}    double* sum = (double*)cv::alignPtr(base_buf,sizeof(double)) + 1;    double* counts = sum + mi + 1;    double** sum_ptr = (double**)(counts + mi);    double L = 0, R = 0, best_val = init_quality, lsum = 0, rsum = 0;    int i, best_subset = -1, subset_i;    for( i = -1; i < mi; i++ )//mi 256,LBP一共256個category        sum[i] = counts[i] = 0;    // calculate sum response and weight of each category of the input var    for( i = 0; i < n; i++ )//對於所有正負樣本    {        int idx = ((cat_labels[i] == 65535) && data->is_buf_16u) ? -1 : cat_labels[i];//對於LBP來說，這個idx都屬於[0,255]        double w = weights[i];        double s = sum[idx] + responses[i]*w;//把weight給乘以(+1)或(-1)求和        double nc = counts[idx] + w;//純粹weight求和        sum[idx] = s;        counts[idx] = nc;    }    // calculate average response in each category    for( i = 0; i < mi; i++ )//對於LBP的256個category    {        R += counts[i];        rsum += sum[i];        sum[i] = fabs(counts[i]) > DBL_EPSILON ? sum[i]/counts[i] : 0;//求category裡的平均值，要是這個category不太重要的話那麼這個值估計會比較小(極端點考慮如果一半是正樣本，一半是負樣本那正好為0)        sum_ptr[i] = sum + i;    }    icvSortDblPtr( sum_ptr, mi, 0 );//把256個值進行升序排序，注意sum_ptr裡存的是sum[i]的地址    // revert back to unnormalized sums    // (there should be a very little loss in accuracy)    for( i = 0; i < mi; i++ )        sum[i] *= counts[i];    for( subset_i = 0; subset_i < mi-1; subset_i++ )    {        int idx = (int)(sum_ptr[subset_i] - sum);//表示排序之後的第subset_i個在排序前是第idx個        double ni = counts[idx];        if( ni > FLT_EPSILON )        {            double s = sum[idx];            lsum += s; L += ni;            rsum -= s; R -= ni;            if( L > FLT_EPSILON && R > FLT_EPSILON )            {                double val = (lsum*lsum*R + rsum*rsum*L)/(L*R);//要賦值給下面的quality，然後外層的每個feature的迴圈裡還要用到                if( best_val < val )                {                    best_val = val;                    best_subset = subset_i;                }            }        }    }    CvDTreeSplit* split = 0;    if( best_subset >= 0 )    {        split = _split ? _split : data->new_split_cat( 0, -1.0f);        split->var_idx = vi;//表示這是利用的第幾個feature        split->quality = (float)best_val;//上層的每個feature的迴圈裡要用來比較好選擇哪個feature        memset( split->subset, 0, (data->max_c_count + 31)/32 * sizeof(int));        for( i = 0; i <= best_subset; i++ )        {            int idx = (int)(sum_ptr[i] - sum);            split->subset[idx >> 5] |= 1 << (idx & 31);//在這裡把256個bin放到8個subset中（8 x 32 = 256），每個subset裡麵包含32位資訊來表示這32個情況是存在與否，//但是我的問題是為啥subset定義的大小是2，根本不夠使的啊。        }    }    return split;}

就是在這裡最後把split->subset[]給寫入資料的。

個人感覺sum[i] = fabs(counts[i]) > DBL_EPSILON ? sum[i]/counts[i] : 0;是[1] 裡的Eq.4，但是double val = (lsum*lsum*R + rsum*rsum*L)/(L*R);這個公式從哪來的不明白，[1]裡對應的應該是Eq.2才對啊。

在xml檔案裡除了subset還有threshold和output值，關於threshold的思路跟haartraing差不多，關於output值（即left->value 和 right->value）是在boost.cpp的calc_node_value(CvDTreeNode *n)中的。主要代碼是：

{        // in case of regression tree:        //  * node value is 1/n*sum_i(Y_i), where Y_i is i-th response,        //    n is the number of samples in the node.        //  * node risk is the sum of squared errors: sum_i((Y_i - <node_value>)^2)        double sum = 0, sum2 = 0, iw;        float* values_buf = (float*)(labels_buf + n);        int* sample_indices_buf = (int*)(values_buf + n);        const float* values = data->get_ord_responses(node, values_buf, sample_indices_buf);        for( i = 0; i < n; i++ )        {            int idx = labels[i];            double w = weights[idx]/*priors[values[i] > 0]*/;double t = values[i];//{-1,1}            rcw[0] += w;            subtree_weights[i] = w;            sum += t*w;            sum2 += t*t*w;        }        iw = 1./rcw[0];        node->value = sum*iw;        node->node_risk = sum2 - (sum*iw)*sum;        // renormalize the risk, as in try_split_node the unweighted formula        // sqrt(risk)/n is used, rather than sqrt(risk)/sum(weights_i)        node->node_risk *= n*iw*n*iw;    }

這部分代碼對應的數學公式是出自哪裡啊？有知道的朋友幫忙告訴一聲。

【1】 Face Detection Based on Multi-Block LBP Representation, L.Zhang