train _cascade 源码阅读之HOG特征

本文讨论OpenCV  train_cascade 级联分类器中的HOG特征实现，HOG特征原理可以参考此文。特征的初始化框架和LBP 特征是一致的，感兴趣可以参考

train_cascade 源码阅读之LBP 特征 中的介绍。 HOG，即Histogram of Oriented Gradient 方向梯度直方图，常用于解决人体目标的检测的图像描述子，用来表达人体，提取人体外形信息和运动信息形成丰富的特征集。
生成过程： 检测窗口--> 归一化图像--> 计算梯度-->对每一个cell块对梯度直方图进行规定权重的投影 --> 对每个重叠block块内的cell进行对比度归一化 --> 把所有block内的直方图向量一起组合成一个大的HOG特征向量。（参考自 blog.sina.com.cn/s/blog_60e6e3d50101bkpn.html）
在HOG特征的操作中，与前述的Haar特征与LBP特征不同的是，初始时矩阵不再是单一的一个sum，而是一个矩阵向量hist，其中含有9个类似sum一样的矩阵，分别存放每一个方向的结果。

[cpp]  view plain  copy

 

1. class CvHOGEvaluator : public CvFeatureEvaluator  
2. {  
3. public:  
4. ……  
5.     virtual void setImage(const cv::Mat& img, uchar clsLabel, int idx);  
6.     virtual float operator()(  
7.             int varIdx,   
8.             int sampleIdx) const;  
9. ……  
10. protected:  
11.     virtual void generateFeatures();  
12.     virtual void integralHistogram(  
13.             const cv::Mat &img,   
14.             std::vector<cv::Mat> &histogram,   
15.             cv::Mat &norm,   
16.             int nbins) const;  
17.     class Feature  
18.     {  
19.     public:  
20.         Feature();  
21.         Feature( int offset, int x, int y, int cellW, int cellH );  
22.         float calc(   
23.                 const std::vector<cv::Mat> &_hists,   
24.                 const cv::Mat &_normSum, size_t y,   
25.                 int featComponent ) const;  
26.         void write( cv::FileStorage &fs ) const;  
27.         void write( cv::FileStorage &fs, int varIdx ) const;  
28.   
29.         cv::Rect rect[N_CELLS]; //cells  
30.   
31.         struct  
32.         {  
33.             int p0, p1, p2, p3;  
34.         } fastRect[N_CELLS];  
35.     };  
36.     std::vector<Feature> features;  
37.   
38.     cv::Mat normSum; //for nomalization calculation (L1 or L2)  
39.     std::vector<cv::Mat> hist;  
40. };  

接下来是初始化积分图操作下标的过程。

[cpp]  view plain  copy

 

1. void CvHOGEvaluator::generateFeatures()  
2. {  
3.     int offset = winSize.width + 1;  
4.     Size blockStep;  
5.     int x, y, t, w, h;  
6.   
7.     for (t = 8; t <= winSize.width/2; t+=8)   
8.         //t = size of a cell. blocksize = 4*cellSize  
9.     {  
10.         blockStep = Size(4,4);  
11.         w = 2*t; //width of a block  
12.         h = 2*t; //height of a block  
13.         for (x = 0; x <= winSize.width - w; x += blockStep.width)  
14.         {  
15.             for (y = 0; y <= winSize.height - h; y += blockStep.height)  
16.             {  
17.                 features.push_back(Feature(offset, x, y, t, t));  
18.             }  
19.         }  
20.         w = 2*t;  
21.         h = 4*t;  
22.         for (x = 0; x <= winSize.width - w; x += blockStep.width)  
23.         {  
24.             for (y = 0; y <= winSize.height - h; y += blockStep.height)  
25.             {  
26.                 features.push_back(Feature(offset, x, y, t, 2*t));  
27.             }  
28.         }  
29.         w = 4*t;  
30.         h = 2*t;  
31.         for (x = 0; x <= winSize.width - w; x += blockStep.width)  
32.         {  
33.             for (y = 0; y <= winSize.height - h; y += blockStep.height)  
34.             {  
35.                 features.push_back(Feature(offset, x, y, 2*t, t));  
36.             }  
37.         }  
38.     }  
39.   
40.     numFeatures = (int)features.size();  
41. }  

t表示cell的尺寸，一个block含有2×2个cell，因此，t需要不大于winSize.width/2，在这里采用了三种不同的形状，block遍历的step是4×4的。输入给Feature的构造参数是偏移量，左上角坐标点，和cell的宽高。
接下来看Feature的构造。

[cpp]  view plain  copy

 

1. CvHOGEvaluator::Feature::Feature( int offset, int x, int y, int cellW, int cellH )  
2. {  
3.     rect[0] = Rect(x, y, cellW, cellH); //cell0  
4.     rect[1] = Rect(x+cellW, y, cellW, cellH); //cell1  
5.     rect[2] = Rect(x, y+cellH, cellW, cellH); //cell2  
6.     rect[3] = Rect(x+cellW, y+cellH, cellW, cellH); //cell3  
7.   
8.     for (int i = 0; i < N_CELLS; i++)  
9.     {  
10.         CV_SUM_OFFSETS(fastRect[i].p0, fastRect[i].p1, fastRect[i].p2, fastRect[i].p3, rect[i], offset);  
11.     }  
12. }  

分别创建了四个cell矩形，CV_SUM_OFFSET宏计算的是矩形上的点在拉成行向量的积分图中的偏移量。
积分图中的坐标算好了，再看积分图的生成过程。

[cpp]  view plain  copy

 

1. void CvHOGEvaluator::  
2. setImage(const Mat &img, uchar clsLabel, int idx)  
3. {  
4.     CV_DbgAssert( !hist.empty());  
5.     CvFeatureEvaluator::setImage( img, clsLabel, idx );  
6.     vector<Mat> integralHist;  
7.     for (int bin = 0; bin < N_BINS; bin++)  
8.     {  
9.         integralHist.push_back(   
10.                     Mat(winSize.height + 1,   
11.                         winSize.width + 1,   
12.                         hist[bin].type(),   
13.                         hist[bin].ptr<float>((int)idx)) );  
14.     }  
15.     Mat integralNorm(  
16.                 winSize.height + 1,   
17.                 winSize.width + 1,   
18.                 normSum.type(),   
19.                 normSum.ptr<float>((int)idx));  
20.     integralHistogram(img, integralHist, integralNorm, (int)N_BINS);  
21. }  

与LBP，Haar相同，新建Mat，传入积分图的数据地址，不同的是这里是矩阵向量，保存9个方向的积分图。这里没有使用OpenCV自带的integral直接计算积分图，而是自行实现了一个。

[cpp]  view plain  copy

 

1. void CvHOGEvaluator::integralHistogram(  
2.         const Mat   &img,   
3.         vector<Mat> &histogram,   
4.         Mat         &norm,   
5.         int         nbins) const  
6. {  
7.     CV_Assert( img.type() == CV_8U || img.type() == CV_8UC3 );  
8.     int x, y, binIdx;  
9.   
10.     Size gradSize(img.size());  
11.     Size histSize(histogram[0].size());  
12.     Mat grad(gradSize, CV_32F);  
13.     Mat qangle(gradSize, CV_8U);  
14.   
15.     AutoBuffer<int> mapbuf(gradSize.width + gradSize.height + 4);  
16.     int* xmap = (int*)mapbuf + 1;  
17.     int* ymap = xmap + gradSize.width + 2;  
18.   
19.     const int borderType = (int)BORDER_REPLICATE;  
20.   
21.     for( x = -1; x < gradSize.width + 1; x++ )  
22.         xmap[x] = borderInterpolate(x, gradSize.width, borderType);  
23.     for( y = -1; y < gradSize.height + 1; y++ )  
24.         ymap[y] = borderInterpolate(y, gradSize.height, borderType);  
25.   
26.     int width = gradSize.width;  
27.     AutoBuffer<float> _dbuf(width*4);  
28.     float* dbuf = _dbuf;  
29.     Mat Dx(1, width, CV_32F, dbuf);  
30.     Mat Dy(1, width, CV_32F, dbuf + width);  
31.     Mat Mag(1, width, CV_32F, dbuf + width*2);  
32.     Mat Angle(1, width, CV_32F, dbuf + width*3);  
33.   
34.     float angleScale = (float)(nbins/CV_PI);  
35.   
36.     for( y = 0; y < gradSize.height; y++ )  
37.     {  
38.         const uchar* currPtr = img.data + img.step*ymap[y];  
39.         const uchar* prevPtr = img.data + img.step*ymap[y-1];  
40.         const uchar* nextPtr = img.data + img.step*ymap[y+1];  
41.         float* gradPtr = (float*)grad.ptr(y);  
42.         uchar* qanglePtr = (uchar*)qangle.ptr(y);  
43.   
44.         for( x = 0; x < width; x++ )  
45.         {  
46.             dbuf[x] = (float)(currPtr[xmap[x+1]] - currPtr[xmap[x-1]]);  
47.             dbuf[width + x] = (float)(nextPtr[xmap[x]] - prevPtr[xmap[x]]);  
48.         }  
49.         cartToPolar( Dx, Dy, Mag, Angle, false );  
50.         for( x = 0; x < width; x++ )  
51.         {  
52.             float mag = dbuf[x+width*2];  
53.             float angle = dbuf[x+width*3];  
54.             angle = angle*angleScale - 0.5f;  
55.             int bidx = cvFloor(angle);  
56.             angle -= bidx;  
57.             if( bidx < 0 )  
58.                 bidx += nbins;  
59.             else if( bidx >= nbins )  
60.                 bidx -= nbins;  
61.   
62.             qanglePtr[x] = (uchar)bidx;  
63.             gradPtr[x] = mag;  
64.         }  
65.     }  
66.     integral(grad, norm, grad.depth());  
67.   
68.     float* histBuf;  
69.     const float* magBuf;  
70.     const uchar* binsBuf;  
71.   
72.     int binsStep = (int)( qangle.step / sizeof(uchar) );  
73.     int histStep = (int)( histogram[0].step / sizeof(float) );  
74.     int magStep = (int)( grad.step / sizeof(float) );  
75.     for( binIdx = 0; binIdx < nbins; binIdx++ )  
76.     {  
77.         histBuf = (float*)histogram[binIdx].data;  
78.         magBuf = (const float*)grad.data;  
79.         binsBuf = (const uchar*)qangle.data;  
80.   
81.         memset( histBuf, 0, histSize.width * sizeof(histBuf[0]) );  
82.         histBuf += histStep + 1;  
83.         for( y = 0; y < qangle.rows; y++ )  
84.         {  
85.             histBuf[-1] = 0.f;  
86.             float strSum = 0.f;  
87.             for( x = 0; x < qangle.cols; x++ )  
88.             {  
89.                 if( binsBuf[x] == binIdx )  
90.                     strSum += magBuf[x];  
91.                 histBuf[x] = histBuf[-histStep + x] + strSum;  
92.             }  
93.             histBuf += histStep;  
94.             binsBuf += binsStep;  
95.             magBuf += magStep;  
96.         }  
97.     }  
98. }  

看完了代码才知道，以前对HOG特征的理解是有偏差的，尤其是在梯度的计算上，犯了严重的想当然的错误。特征计算完成后，调用integral计算平方积分图，再根据角度，将幅值放到每个积分直方图中。
最后通过如下方式调用计算HOG特征，并进行归一化。

[cpp]  view plain  copy

 

1. inline float CvHOGEvaluator::Feature::calc(   
2.         const std::vector<cv::Mat>& _hists,   
3.         const cv::Mat& _normSum, size_t y,   
4.         int featComponent ) const  
5. {  
6.     float normFactor;  
7.     float res;  
8.   
9.     int binIdx = featComponent % N_BINS;  
10.     int cellIdx = featComponent / N_BINS;  
11.   
12.     const float *phist = _hists[binIdx].ptr<float>((int)y);  
13.     res = phist[fastRect[cellIdx].p0]   
14.             - phist[fastRect[cellIdx].p1]   
15.             - phist[fastRect[cellIdx].p2]   
16.             + phist[fastRect[cellIdx].p3];  
17.   
18.     const float *pnormSum = _normSum.ptr<float>((int)y);  
19.     normFactor = (float)(pnormSum[fastRect[0].p0]   
20.             - pnormSum[fastRect[1].p1]   
21.             - pnormSum[fastRect[2].p2]   
22.             + pnormSum[fastRect[3].p3]);  
23.     res = (res > 0.001f) ? ( res / (normFactor + 0.001f) ) : 0.f;   
24.     //for cutting negative values, which apper due to floating precision  
25.   
26.     return res;  
27. }