research►note►computer vision

HOG不正传

闲话历史

HOG(histogram of oriented gradients)是object detection中非常popular的特征，对于行人检测尤其有效，能够较好地捕捉行人轮廓。最早由 Navneet Dalal和Bill Triggs在2005年的CVPR论文Histograms of Oriented Gradients for Human Detection中提出（这篇神文的应用量已经达到6000+次了）。

所谓“长江后浪推前浪，一代更比一代浪(强)！”. Brown University的pff使用HOG特征加上DPM模型和latent SVM学习，多次在VOC challenge中夺魁 (btw, DPM是他在University of Chicago的工作) ，最终获得lifetime achievement的title(新闻见这里).

又所谓“不开源代码的工作都是耍流氓”，作为一代神作，pff和他的合作者们孜孜不倦地迭代了5个版本的DPM software package，引用另一个我很欣赏的计算机视觉研究者Tomasz Malisiewicz (顺便推荐一下此人的blog, 有很多干货)的话来说，”A key ingredient to successful object recognition research is a powerful codebase, which you will hopefully one day outgrow and/or extend. “（原文见Tomasz的日志）

DPM software package如此流行，所以本文结合DPM PAMI论文中的理论及软件包中features.cc的代码介绍HOG特征的计算。关于OpenCV的HOG源码分析也有不少资源，供使用OpenCV的研究者参考。

Why is HOG effective?

为啥HOG特征这么有效呢？ 手工设计一个特征，比如DOG, HOG, LBP等等（感兴趣的童鞋可以戳Zhu Songchun主页上关于CV研究的漫画），在我看来HOG特征的设计如同飞来神笔，仍然是一件非常神秘的事情。对着这个问题的解释主要参考了HOG的wiki页面

HOG的基本思想是，图像中局部物体的轮廓和形状可以用灰度值梯度强度和方向的分布描述。将图像划分成小块的连通区域(称作cell, 这是一种“sliding window”的方法), 计算cell内像素点的梯度方向的直方图(histogram of gradient directions). 为了使特征对光照和阴影保持不变性，将几个cell聚合成一个block, 对直方图进行contrast-normalization.
相对于其他特征，HOG有以下两点优势

  • 对局部区域进行操作，对geometric and photometric transformation保持不变性
  • 对于行人，只要目标保持直立姿态(upright), HOG特征不受身体运动的影响，因此HOG特征对行人检测效果特别好

题外话 关于如何学出一个特征这个一度像黑匣子一样神秘的问题，Deep Learning研究在该问题上取得了一定进展，不过我对此方向不太了解，且这不是本文的话题，保持持续关注，暂且按下不表。pluskid大神关于deep learning发表了一篇blog, 不明觉厉，值得一看。

灌了好些水，也该上点干货了，请猛烈点击下方的Read on按钮继续阅读 XD.

Pixel-level feature maps

令$\theta(x, y)$和$\gamma(x, y)$分别表示图像中$(x, y)$坐标处的方向(orientation)和强度(magnitude). 使用模板[-1, 0, 1]及其转置计算梯度。对于彩色图像，计算每个通道的梯度，选择其中的最大梯度值作为该点的梯度值。

将gradient orientation离散化到$p$个值，分别记contrast sensitive和insensitive为$$B_1$$和$$B_2$$，简记为$$B$$. 其中

$$
\begin{equation}
\begin{cases}
B_1(x, y) = round(\frac{p\theta(x, y)}{2 \pi})\ mod\ p \
B_2(x, y) = round(\frac{p\theta(x, y)}{\pi})\ mod\ p
\end{cases}
\end{equation}
$$

对每个像素点$(x, y)$定义feature map, 其中$b \in {0, \ldots, p-1}$

$$
\begin{equation}
F(x, y) =
\begin{cases}
r(x, y) & if\ b = B(x, y) \
0 & otherwise
\end{cases}
\end{equation}
$$

注 Feature map是一个p维稀疏向量，除b外的值都是0. 将feature map看作oriented edge map, 根据gradient orientation选择p个通道中的一个，

Spatial Aggregation

记$F$为大小为$w \times h$的图像的feature map, 将pixel-level feature map $F$聚合成cell-based feature map $C$, 每个cell是一个$k \times k$的正方形. 空间聚合能够使特征对微小形变具有不变性，同时减小feature map的维度。
聚合方法采用”soft binning”的方法，用双线性插值法，将每个pixel feature map分配到与其相邻的4个cell中。具体方法见下节

Normalization and truncation

  • Normalization factor的计算

$$
\begin{equation}
N_{\delta, \gamma} = (||C(i, j)||^2 + ||C(i+\delta, j)||^2 + ||C(i, j+\gamma)||^2 + ||C(i+\delta, j+\gamma)||^2)^{\frac{1}{2}}
\end{equation}
$$

其中$\delta, \gamma \in {-1, 1}$
每个normalization factor计算包括$C(i, j)$的4个cell领域构成的block的gradient energy, $N_{-1, -1}$的block构成如下图所示，绿色区域为$C(i, j)$, 其他blockd的构成可依此类推。



  • 定义Trucation函数$T_{\alpha}(v)$, 将向量$v$中的任一元素$v(i)$与$\alpha$比较，取其中的较小值。

最终HoG feature map将经过truncation和normalization操作的的cell-based feature map $C$连接在一起。

$$
\begin{equation}
H(i, j) =
\begin{pmatrix}
T{\alpha}(C(i, j)/N{-1, -1}(i, j)) \
T{\alpha}(C(i, j)/N{+1, -1}(i, j)) \
T{\alpha}(C(i, j)/N{+1, +1}(i, j)) \
T{\alpha}(C(i, j)/N{-1, +1}(i, j))
\end{pmatrix}
\end{equation}
$$

常用的参数为cell size $k = 8$, $\alpha = 0.2$, $p = 9$.

代码分析

features.cc是DPM software package中用于计算HoG特征的c++代码，package主要用MATLAB开发，计算量较大的部分如HoG特征和卷积的计算使用c++语言开发。

  • mex函数入口
    函数原型为F = features(image, bin)，第一个输入参数为三通道彩色图像，灰度值必须为double型(MATLAB的默认精度)，第二个输入参数为cell size的边长，即上文中的$k$, 返回值为输入图像image的HoG特征。
1
2
3
4
5
6
7
8
9
10
// matlab entry point
// F = features(image, bin)
// image should be color with double values
void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) { 
  if (nrhs != 2)
    mexErrMsgTxt("Wrong number of inputs"); 
  if (nlhs != 1)
    mexErrMsgTxt("Wrong number of outputs");
  plhs[0] = process(prhs[0], prhs[1]);
}
  • HoG特征计算主函数process
    cache some values for looking up. 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
// small value, used to avoid division by zero
#define eps 0.0001

// unit vectors used to compute gradient orientation
double uu[9] = {1.0000, 
		0.9397, 
		0.7660, 
		0.500, 
		0.1736, 
		-0.1736, 
		-0.5000, 
		-0.7660, 
		-0.9397};
double vv[9] = {0.0000, 
		0.3420, 
		0.6428, 
		0.8660, 
		0.9848, 
		0.9848, 
		0.8660, 
		0.6428, 
		0.3420};

函数原型mxArray *process(const mxArray *mximage, const mxArray *mxsbin), mxArray为mex内建类型，mximage为输入图像，mxsbin为cell size边长

1
2
3
4
// main function:
// takes a double color image and a bin size 
// returns HOG features
mxArray *process(const mxArray *mximage, const mxArray *mxsbin)

将mex内建类型转换成c++数据类型

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
double *im = (double *)mxGetPr(mximage); // im指针指向输入图像
const int *dims = mxGetDimensions(mximage); // 输入图像的维度信息[h, w, 深度]
// 图像必须是三通道彩色图像，数据精度为double
if (mxGetNumberOfDimensions(mximage) != 3 ||
    dims[2] != 3 ||
    mxGetClassID(mximage) != mxDOUBLE_CLASS)
  mexErrMsgTxt("Invalid input");

int sbin = (int)mxGetScalar(mxsbin); // cell size的边长

// memory for caching orientation histograms & their norms
int blocks[2];
blocks[0] = (int)round((double)dims[0]/(double)sbin); // cell的个数 h/k
blocks[1] = (int)round((double)dims[1]/(double)sbin); // cell个数 w/k
// 为计算pixel feature map和normalized feature map开辟空间
// p = 9, contrast sensitive B1和contrast insensitive B2，共开辟空间18
double *hist = (double *)mxCalloc(blocks[0]*blocks[1]*18, sizeof(double));
double *norm = (double *)mxCalloc(blocks[0]*blocks[1], sizeof(double));

// memory for HOG features
int out[3];
// out[0], out[1]减2的原因：上下边缘无法计算4领域normalization
out[0] = max(blocks[0]-2, 0);	//输出的y维度
out[1] = max(blocks[1]-2, 0);	//输出的x维度
out[2] = 27+4+1;  //每个block的特征维度
mxArray *mxfeat = mxCreateNumericArray(3, out, mxDOUBLE_CLASS, mxREAL);
double *feat = (double *)mxGetPr(mxfeat);

int visible[2];
visible[0] = blocks[0]*sbin;	
visible[1] = blocks[1]*sbin;

计算每一点$(x, y)$的pixel-level feature map 

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
 // 遍历每个visible pixel (x, y)
 for (int x = 1; x < visible[1]-1; x++) {
   for (int y = 1; y < visible[0]-1; y++) {
     // 计算每个色彩通道的gradient orientation
     // first color channel
     double *s = im + min(x, dims[1]-2)*dims[0] + min(y, dims[0]-2);
     double dy = *(s+1) - *(s-1);
     double dx = *(s+dims[0]) - *(s-dims[0]);
     double v = dx*dx + dy*dy;

     // second color channel
     s += dims[0]*dims[1];
     double dy2 = *(s+1) - *(s-1);
     double dx2 = *(s+dims[0]) - *(s-dims[0]);
     double v2 = dx2*dx2 + dy2*dy2;

     // third color channel
     s += dims[0]*dims[1];
     double dy3 = *(s+1) - *(s-1);
     double dx3 = *(s+dims[0]) - *(s-dims[0]);
     double v3 = dx3*dx3 + dy3*dy3;

     // pick channel with strongest gradient
     if (v2 > v) {  // 把最强的边缘保持在v，dx，dy里
v = v2;
dx = dx2;
dy = dy2;
     } 
     if (v3 > v) {
v = v3;
dx = dx3;
dy = dy3;
     }

     // snap to one of 18 orientations
     double best_dot = 0;
     int best_o = 0;
     for (int o = 0; o < 9; o++) {
	double dot = uu[o]*dx + vv[o]*dy; // 在9个方向上，正负就是18个方向
	if (dot > best_dot) {		//在那个方向上的响应最大，结果保存在best_dot里
	  best_dot = dot;
	  best_o = o;
	} else if (-dot > best_dot) {
	  best_dot = -dot;
	  best_o = o+9;
	}
     }
     
     // add to 4 histograms around pixel using linear interpolation
     double xp = ((double)x+0.5)/(double)sbin - 0.5;
     double yp = ((double)y+0.5)/(double)sbin - 0.5;
     int ixp = (int)floor(xp);
     int iyp = (int)floor(yp);
     double vx0 = xp-ixp;  //在block里的位置比例
     double vy0 = yp-iyp;  
     double vx1 = 1.0-vx0;  //在block里另外一边的位置比例
     double vy1 = 1.0-vy0;
     v = sqrt(v);

     if (ixp >= 0 && iyp >= 0) { //像素对应那个block 一个像素要投在4个block里
*(hist + ixp*blocks[0] + iyp + best_o*blocks[0]*blocks[1]) += 
  vx1*vy1*v;	//在block里的位置和强度乘积
     }

     if (ixp+1 < blocks[1] && iyp >= 0) {
*(hist + (ixp+1)*blocks[0] + iyp + best_o*blocks[0]*blocks[1]) += 
  vx0*vy1*v;
     }

     if (ixp >= 0 && iyp+1 < blocks[0]) {
*(hist + ixp*blocks[0] + (iyp+1) + best_o*blocks[0]*blocks[1]) += 
  vx1*vy0*v;
     }

     if (ixp+1 < blocks[1] && iyp+1 < blocks[0]) {
*(hist + (ixp+1)*blocks[0] + (iyp+1) + best_o*blocks[0]*blocks[1]) += 
  vx0*vy0*v;
     }
   }
 }
   // compute energy in each block by summing over orientations
   // 计算normalization factors
 for (int o = 0; o < 9; o++) { //每个方向分别做，0和9 1和10， 每次处理2个方向
   double *src1 = hist + o*blocks[0]*blocks[1];
   double *src2 = hist + (o+9)*blocks[0]*blocks[1];
   double *dst = norm;
   double *end = norm + blocks[1]*blocks[0];
   while (dst < end) {
     *(dst++) += (*src1 + *src2) * (*src1 + *src2); //成对方向相加后平方，累计所有的9对方向
     src1++;
     src2++;
   }
 }

计算HoG特征, concatenating normalized cell-based feature map

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
// compute features
  for (int x = 0; x < out[1]; x++) {
    for (int y = 0; y < out[0]; y++) {
      double *dst = feat + x*out[0] + y;      
      double *src, *p, n1, n2, n3, n4; // n1-n4为4个normalization factors
		//统计 16个block的综合，归一化，每4个一组
      p = norm + (x+1)*blocks[0] + y+1; //y+1， x+1block里的norm ， 
      n1 = 1.0 / sqrt(*p + *(p+1) + *(p+blocks[0]) + *(p+blocks[0]+1) + eps); //4个block里的和
      p = norm + (x+1)*blocks[0] + y;
      n2 = 1.0 / sqrt(*p + *(p+1) + *(p+blocks[0]) + *(p+blocks[0]+1) + eps);
      p = norm + x*blocks[0] + y+1;
      n3 = 1.0 / sqrt(*p + *(p+1) + *(p+blocks[0]) + *(p+blocks[0]+1) + eps);
      p = norm + x*blocks[0] + y;      
      n4 = 1.0 / sqrt(*p + *(p+1) + *(p+blocks[0]) + *(p+blocks[0]+1) + eps);

      double t1 = 0;
      double t2 = 0;
      double t3 = 0;
      double t4 = 0;

      // contrast-sensitive features	
      src = hist + (x+1)*blocks[0] + (y+1);
      for (int o = 0; o < 18; o++) {
      // truncation操作, h1-h4为H(i, j)的4个分量
	double h1 = min(*src * n1, 0.2);
	double h2 = min(*src * n2, 0.2);
	double h3 = min(*src * n3, 0.2);
	double h4 = min(*src * n4, 0.2);
	*dst = 0.5 * (h1 + h2 + h3 + h4);
	// t1-t4 cache for computing texture
	t1 += h1;
	t2 += h2;
	t3 += h3;
	t4 += h4;
	dst += out[0]*out[1];
	src += blocks[0]*blocks[1];
      }

      // contrast-insensitive features
      src = hist + (x+1)*blocks[0] + (y+1);
      for (int o = 0; o < 9; o++) {
        double sum = *src + *(src + 9*blocks[0]*blocks[1]);
        double h1 = min(sum * n1, 0.2);
        double h2 = min(sum * n2, 0.2);
        double h3 = min(sum * n3, 0.2);
        double h4 = min(sum * n4, 0.2);
        *dst = 0.5 * (h1 + h2 + h3 + h4);
        dst += out[0]*out[1];
        src += blocks[0]*blocks[1];
      }

      // texture features
      *dst = 0.2357 * t1;
      dst += out[0]*out[1];
      *dst = 0.2357 * t2;
      dst += out[0]*out[1];
      *dst = 0.2357 * t3;
      dst += out[0]*out[1];
      *dst = 0.2357 * t4;

      // truncation feature
      dst += out[0]*out[1];
      *dst = 0;
    }
  }

  // 释放空间，返回特征计算结果
  mxFree(hist);
  mxFree(norm);
  return mxfeat;