本文深入解析HOG特征的原理及提取流程，对比SIFT特征，介绍HOG在图像局部梯度方向信息统计上的优势，以及在行人检测和目标跟踪领域的应用。
HOG 特征 — 学习理解
Email: animatex.deng@gmail.com
Reference:
- Paper: [1] Histograms of Oriented Gradients for Human Detection
- Paper: [2] Finding People in Images and Videos
- Blog: [3] https://blog.youkuaiyun.com/zhazhiqiang/article/details/21047207
- Blog: [4] https://blog.youkuaiyun.com/langb2014/article/details/45674725
- Zhihu: [5] https://zhuanlan.zhihu.com/p/33059421
0 HOG特征 对比 SIFT特征
HOG - Histograms of Oriented Gradients; SIFT - Scale-invariant feature transform
最近看跟踪相关论文时，对handcraft特征很感兴趣，那就顺势理解一波。先看一两张表来看看HOG特征的优劣和SIFT特征对比。需要注意的是HOG特征不是SIFT特征中的子类。
| Feature | 优点 | 缺点 | 特征适用范围 |
|---|---|---|---|
| SIFT | ① 对于旋转、尺度缩放、亮度变化保持不变性；② 抗遮挡 | 计算量大 | ① 图像匹配；② 三维建模 |
| HOG | ① 减弱了光照颜色对图像造成的影响，使得图像所需要表征的数据的维度降低；② 描述局部形状信息；③ 一定程度抑制平移和旋转带来的影响。 | ① 描述子生成过程冗长，速度较慢，实时性较差；② 不抗遮挡；③ 对噪声及其敏感 | ① 行人检测；② 目标跟踪 |
| LBP | ① 对光照不敏感；② 运算速度快 | 对方向信息敏感 | ① 人脸识别；② 照片分类 |
HOG和SIFT都是描述子，以及由于在具体操作上有很多相似的步骤，所以致使很多人误认为HOG是SIFT的一种，其实两者在使用目的和具体处理细节上是有很大的区别的。HOG与SIFT的主要区别如下：
- SIFT是基于关键点特征向量的描述。
- HOG是将图像均匀的分成相邻的小块，然后在所有的小块内统计梯度直方图。
- SIFT需要对图像尺度空间下对像素求极值点，而HOG中不需要。
- SIFT一般有两大步骤，第一个步骤对图像提取特征点，而HOG不会对图像提取特征点。
1 HOG特征本质
Locally normalised histogram of gradient orientation in dense overlapping grids.
Histogram of Oriented Gradient descriptors provide a dense overlapping description of image regions.
统计图像局部梯度方向信息（直方图）作为局部区域的表征。
2 HOG特征提取流程
2.1 基本信息
如图所示，灰色为输入的图像，蓝色为检测窗口大小，红色为block大小，绿色为cell大小。
win_size: 64 * 128
block_size: 16 * 16
cell_size: 8 * 8
hist_bin: 9 ---> [0 : 20 : pi] 每一个cell统计9个方向的梯度直方图，构成'描述子'向量长度为9
slide stride: cell_size
2.2 算法流程
① 图像灰度化
颜色信息对HOG特征没有贡献，输入彩色图像转换为灰度图。
② 标准化gamma空间
为了减少光照影响，实测貌似没有太大的用处，实际使用中在这一步后会对进来的图片resize到window的Size，参考后续代码。
I
(
x
,
y
)
=
I
(
x
,
y
)
g
a
m
m
a
\mathrm{I(x,y) = I(x,y)^{gamma}}
I(x,y)=I(x,y)gamma
③ 计算梯度大小和方向
Operator，x方向：[-1 0 1]； y方向：[-1 0 1]'
G
(
x
,
y
)
=
G
x
2
+
G
y
2
\mathrm{G(x,y) = \sqrt{G_x^2 + G_y^2}}
G(x,y)=Gx2+Gy2
α ( x , y ) = t a n − 1 ( G y G x ) \mathrm{\alpha(x,y) = tan^{-1}(\frac{G_y}{G_x})} α(x,y)=tan−1(GxGy)
④ 将图像分割成多个cell
由前述可知，将图像划分成多个cell（8 * 8 pixel）。
⑤ 生成每个cell的梯度方向直方图 [重点]
统计局部图像梯度信息并进行量化（或称为编码），得到局部图像区域的特征描述向量。
前面已经将图像划分成多个cell（8 * 8），这里默认使用R-HOG（矩形HOG），且采用9个bin的直方图来统计每个cell的64个像素的梯度信息。9个bin实际上就是用9个块均分360°。如下图所示。
例如：如果这个像素的梯度方向是20-40度，直方图第2个bin的计数就加一，这样，对cell内每个像素用梯度方向在直方图中进行加权投影（映射到固定的角度范围），就可以得到这个cell的梯度方向直方图了，就是该cell对应的9维特征向量（因为有9个bin）。梯度大小就是作为投影的权值的。例如说：这个像素的梯度方向是20-40度，然后它的梯度大小是2（假设啊），那么直方图第2个bin的计数就不是加一了，而是加二（假设啊）。
单元格Cell中的每一个像素点都为某个基于方向的直方图通道（orientation-based histogram channel）投票。投票是采取加权投票（weighted voting）的方式，即每一票都是带权值的，这个权值是根据该像素点的梯度幅度计算出来。可以采用幅值本身或者它的函数来表示这个权值，实际测试表明： 使用幅值来表示权值能获得最佳的效果，当然，也可以选择幅值的函数来表示，比如幅值的平方根（square root）、幅值的平方（square of the gradient magnitude）、幅值的截断形式（clipped version of the magnitude）等。根据Dalal等人论文的测试结果，采用梯度幅值量级本身得到的检测效果最佳，使用量级的平方根会轻微降低检测结果，而使用二值的边缘权值表示会严重降低效果。
通过下面几幅图来增强理解。
左上图中梯度方向矩阵中可以看到角度是0 - 180度，不是0 - 360度，这种被称之为"无符号"梯度
("unsigned" gradients)。因为一个梯度和它的负数是用同一个数字表示的，也就是说一个梯度的箭头以及它旋转180度之后的箭头方向被认为是一样的。那为什么不用0 - 360度的表示呢？在事件中发现unsigned gradients比signed gradients在行人检测任务中效果更好。一些HOG的实现中可以让你指定signed gradients。
来看看这个过程，使用9个bin来存放统计结果，其中左上图为梯度方向计算结果，右上图为梯度大小图。先看蓝色圈里的梯度方向为80°。落在[70 90]区间内。对应的权重为梯度大小2，因此在第5个bin中加2。同理在红圈中的10°刚好落在20和40之间，这里将第一和第二bin权重各占2。
这里有个细节要注意，如果一个角度大于160度，也就是在160-180度之间，我们知道这里角度0，180度是一样的，所以在下面这个例子里，像素的角度为165度的时候，由双线性插值简单计算，weight1 = (165-160)/20 = 0.25, weight2 = (180 - 165) / 20 = 0.75，权重分别为0.25 * 85 --> bin1, 0.75 * 85 --> bin9。把幅值按照比例放到0和160的bin里面去，如下图所示。
结果如下：
**实际是怎么做的呢？**加权采用三线性插值方法，即将当前像素的梯度方向大小、像素在cell中的x坐标与y坐标这三个值来作为插值权重，而被用来插入的值为像素的梯度幅值。采用三线性插值的好处在于：避免了梯度方向直方图在cell边界和梯度方向量化的bin边界处的突然变化。这里把链接中的关键点拿出来。参考paper2第131页。
为什么要用插值呢？ 区域混叠造成的，混叠效应会在特征向量中产生突变。在这种情况下就需要采用插值算法对计算结果进行修正。在HOG特征提取方法中，位于不同cell 交界处的像素如果只对所在的cell 进行投影同样会对其他区域产生混叠效应，此时需要采用三维线性插值的方法对梯度向量直方图进行修正。插值运算在图像放大中应用较多，这里借用插值思想对累积直方图进行修正，即将某个像素点的梯度幅值以不同的权重累加到相应的bin上。确切地说，这是一种权值分配，但多数文献中仍然延续“插值”这一说法。由于提取HOG特征时需要在两个位置坐标(x, y)和一个方向坐标( θ ) 共三个维度上进行插值运算，因此称为三线性插值。
首先，以二维情况为例说明线性插值的思想。在两个维度上进行线性插值称为“双线性插值”。 双线性插值在图像放大或图像旋转中应用广泛，它是利用周围4 个邻点的灰度值在两个方向上作线性内插以得到待采样点的灰度值，即根据待采样点与相邻点的距离确定相应的权值计算出待采样点的灰度值。如下图所示，点P为待采样点，Q11、Q12、Q21、Q22为点P 的四个相邻点，用线性插值法对P 点进行插值运算的数学表达式为：
f ( R 1 ) ≈ x 2 − x x 2 − x 1 f ( Q 11 ) + x − x 1 x 2 − x 1 f ( Q 21 ) ⟵ R 1 = ( x , y 1 ) \mathrm{f(R_1)\approx \frac{x_2-x}{x_2-x_1}f(Q_{11}) + \frac{x-x_1}{x_2-x_1}f(Q_{21}) \longleftarrow R1 = (x,y_1)} f(R1)≈x2−x1x2−xf(Q11)+x2−x1x−x1f(Q21)⟵R1=(x,y1)
f ( R 2 ) ≈ x 2 − x x 2 − x 1 f ( Q 12 ) + x − x 1 x 2 − x 1 f ( Q 22 ) ⟵ R 2 = ( x , y 2 ) \mathrm{f(R_2)\approx \frac{x_2-x}{x_2-x_1}f(Q_{12}) + \frac{x-x_1}{x_2-x_1}f(Q_{22}) \longleftarrow R2 = (x,y_2)} f(R2)≈x2−x1x2−xf(Q12)+x2−x1x−x1f(Q22)⟵R2=(x,y2)
f ( P ) ≈ y 2 − y y 2 − y 1 f ( R 1 ) + y − y 1 y 2 − y 1 f ( R 2 ) \mathrm{f(P)\approx \frac{y_2-y}{y_2-y_1}f(R_1) + \frac{y-y_1}{y_2-y_1}f(R_{2})} f(P)≈y2−y1y2−yf(R1)+y2−y1y−y1f(R2)
根据上述插值思想，可在各像素的梯度方向上进行加权运算，如下图所示，将区间[ 0°, 180°] 以20°为一个区间划分，每个小区间以中心角度作为直方图的中心数值，假设要对梯度方向为15°的像素点进行处理，显然15与以10和30为中心的直方图最近，应该将该点的梯度幅值加权累加到这两个直方图上，权重分别为(30 - 15) / 20 = 0. 75 和(15 - 10) / 20 = 0. 25。
归纳为公式：
{
h
(
x
1
)
⟵
h
(
x
1
)
+
ω
(
1
−
x
−
x
1
b
)
h
(
x
2
)
⟵
h
(
x
2
)
+
ω
(
x
−
x
1
b
)
\left\{ \begin{aligned} \mathrm{h(x_1) \longleftarrow h(x_1) + \omega(1-\frac{x-x_1}{b})} \\ \mathrm{h(x_2) \longleftarrow h(x_2) + \omega(\frac{x-x_1}{b})} \\ \end{aligned} \right.
⎩⎪⎨⎪⎧h(x1)⟵h(x1)+ω(1−bx−x1)h(x2)⟵h(x2)+ω(bx−x1)
同理运用线性插值方法在各个像素的位置上进行加权运算，如下图所示，左图中的方框处为待处理像素点，它位于block中的C0 单元中，利用该点与四个cell 中的中心像素点（图中4 个圆点）的距离计算权值，将待处理像素点的梯度幅值分别加权累加到C0、C1、C2、C3 中相应的直方图上。
综合考虑，在两个位置坐标(x, y) 和一个方向坐标( θ )上进行三线性插值，关键要解决的问题是应该在哪些 bin上进行加权累加，累加时权值又是多少。将一个像素点处的梯度幅值加权分配到4 个cell 中与该点梯度方向相邻的2 个bin上。看一下paper中怎么说。
Let w at the 3-D point x = [x, y, z] be the weight to be interpolated.Let x1 and x2 be the two corner vectors of the histogram cube containing x, where in each component x1 <= x <= x2. Assume that the bandwidth of the histogram along the x, y and z axis is given by b = [bx, by, bz]. Trilinear interpolation distributes the weight w to the 8 surroundingbin centres as follows.
点x的x, y, z是要插值的权重，令x1和x2是包含x的直方图立方体的两个角矢量（空间位置），其中在每个分量中x1 <= x <= x2。b为直方图的间距，bx = by = 8, bz默认为20。三线性插值将权重w分配给8个周围的箱中心。按照下面的公式修正直方图向量：
{
h
(
x
1
,
y
1
,
z
1
)
⟵
h
(
x
1
,
y
1
,
z
1
)
+
ω
(
1
−
x
−
x
1
b
x
)
(
1
−
y
−
y
1
b
y
)
(
1
−
z
−
z
1
b
z
)
h
(
x
1
,
y
1
,
z
2
)
⟵
h
(
x
1
,
y
1
,
z
2
)
+
ω
(
1
−
x
−
x
1
b
x
)
(
1
−
y
−
y
1
b
y
)
(
1
−
z
−
z
2
b
z
)
h
(
x
1
,
y
2
,
z
1
)
⟵
h
(
x
1
,
y
2
,
z
1
)
+
ω
(
1
−
x
−
x
1
b
x
)
(
y
−
y
1
b
y
)
(
1
−
z
−
z
1
b
z
)
h
(
x
2
,
y
1
,
z
1
)
⟵
h
(
x
2
,
y
1
,
z
1
)
+
ω
(
x
−
x
1
b
x
)
(
1
−
y
−
y
1
b
y
)
(
1
−
z
−
z
1
b
z
)
h
(
x
1
,
y
2
,
z
2
)
⟵
h
(
x
1
,
y
2
,
z
2
)
+
ω
(
1
−
x
−
x
1
b
x
)
(
y
−
y
1
b
y
)
(
z
−
z
1
b
z
)
h
(
x
2
,
y
1
,
z
2
)
⟵
h
(
x
2
,
y
1
,
z
2
)
+
ω
(
x
−
x
1
b
x
)
(
1
−
y
−
y
1
b
y
)
(
z
−
z
1
b
z
)
h
(
x
2
,
y
2
,
z
1
)
⟵
h
(
x
2
,
y
2
,
z
1
)
+
ω
(
x
−
x
1
b
x
)
(
y
−
y
1
b
y
)
(
1
−
z
−
z
1
b
z
)
h
(
x
2
,
y
2
,
z
2
)
⟵
h
(
x
2
,
y
2
,
z
2
)
+
ω
(
x
−
x
1
b
x
)
(
y
−
y
1
b
y
)
(
z
−
z
1
b
z
)
\left\{ \begin{aligned} \mathrm{h(x_1,y_1,z_1) \longleftarrow h(x_1,y_1,z_1) + \omega(1-\frac{x-x_1}{b_x})(1-\frac{y-y_1}{b_y})(1-\frac{z-z_1}{b_z})} \\ \mathrm{h(x_1,y_1,z_2) \longleftarrow h(x_1,y_1,z_2) + \omega(1-\frac{x-x_1}{b_x})(1-\frac{y-y_1}{b_y})(1-\frac{z-z_2}{b_z})} \\ \mathrm{h(x_1,y_2,z_1) \longleftarrow h(x_1,y_2,z_1) + \omega(1-\frac{x-x_1}{b_x})(\frac{y-y_1}{b_y})(1-\frac{z-z_1}{b_z})} \\ \mathrm{h(x_2,y_1,z_1) \longleftarrow h(x_2,y_1,z_1) + \omega(\frac{x-x_1}{b_x})(1-\frac{y-y_1}{b_y})(1-\frac{z-z_1}{b_z})} \\ \mathrm{h(x_1,y_2,z_2) \longleftarrow h(x_1,y_2,z_2) + \omega(1-\frac{x-x_1}{b_x})(\frac{y-y_1}{b_y})(\frac{z-z_1}{b_z})} \\ \mathrm{h(x_2,y_1,z_2) \longleftarrow h(x_2,y_1,z_2) + \omega(\frac{x-x_1}{b_x})(1-\frac{y-y_1}{b_y})(\frac{z-z_1}{b_z})} \\ \mathrm{h(x_2,y_2,z_1) \longleftarrow h(x_2,y_2,z_1) + \omega(\frac{x-x_1}{b_x})(\frac{y-y_1}{b_y})(1-\frac{z-z_1}{b_z})} \\ \mathrm{h(x_2,y_2,z_2) \longleftarrow h(x_2,y_2,z_2) + \omega(\frac{x-x_1}{b_x})(\frac{y-y_1}{b_y})(\frac{z-z_1}{b_z})} \\ \end{aligned} \right.
⎩⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎧h(x1,y1,z1)⟵h(x1,y1,z1)+ω(1−bxx−x1)(1−byy−y1)(1−bzz−z1)h(x1,y1,z2)⟵h(x1,y1,z2)+ω(1−bxx−x1)(1−byy−y1)(1−bzz−z2)h(x1,y2,z1)⟵h(x1,y2,z1)+ω(1−bxx−x1)(byy−y1)(1−bzz−z1)h(x2,y1,z1)⟵h(x2,y1,z1)+ω(bxx−x1)(1−byy−y1)(1−bzz−z1)h(x1,y2,z2)⟵h(x1,y2,z2)+ω(1−bxx−x1)(byy−y1)(bzz−z1)h(x2,y1,z2)⟵h(x2,y1,z2)+ω(bxx−x1)(1−byy−y1)(bzz−z1)h(x2,y2,z1)⟵h(x2,y2,z1)+ω(bxx−x1)(byy−y1)(1−bzz−z1)h(x2,y2,z2)⟵h(x2,y2,z2)+ω(bxx−x1)(byy−y1)(bzz−z1)
接上，x、y 轴表征像素点的空间位置，z 轴表征该点的梯度方向（即θ ）。对于待处理像素点(x, y) ，设其梯度幅值为ω ，梯度方向为z，z 1 和z 2 分别是与之最邻近的两个bin的中点坐标。下图所示为三线性插值计算梯度方向直方图向量的示意图：
左图中的方框处为待处理像素点，计算 block的每个 cell 中与该点梯度方向相邻的 2 个bin，共计 8 个直方图柱上的权值，将该点的梯度幅值进行加权累加，即形成 block中的梯度方向直方图。
⑥ 将cell组合成block，并归一化block梯度直方图 [重点]
前述可知，HOG提取的是局部特征，由于局部光照的变化以及前景-背景对比度的变化，使得梯度强度的变化范围非常大。这就需要对梯度强度做归一化。归一化能够进一步地对光照、阴影和边缘进行压缩。
如何归一化？
(6-1) 将多个临近的cell组合成一个block块，然后求其梯度方向直方图向量；
(6-2) 采用L2-Norm with Hysteresis threshold方式进行归一化，即将直方图向量中bin值的最大值限制为0.2以下，然后再重新归一化一次；
**注意：block之间的是“共享”的，也即是说，一个cell会被多个block“共享”。**另外，每个“cell”在被归一化时都是“block”independent的，也就是说每个cell在其所属的block中都会被归一化一次，得到一个vector。这就意味着：每一个单元格的特征会以不同的结果多次出现在最后的特征向量中。
(6-3) 四种归一化的方法：
采用了四中不同的方法对区间进行归一化，并对结果进行了比较。引入v表示一个还没有被归一化的向量，它包含了给定区间（block）的所有直方图信息。||Vk||表示v的k阶范数，这里的k去1、2。用e (matlab 中可用eps代替) 表示一个很小的常数。 (L2 norm就是平方和开方，即欧氏距离；L1 norm就是绝对值相加，又称曼哈顿距离)。这时，归一化因子可以表示如下：
- L2-norm, v ⟵ v / ∥ v ∥ 2 2 + ϵ 2 \mathrm{v \longleftarrow v / \sqrt{\|v\|_{2}^{2} + \epsilon^2}} v⟵v/∥v∥22+ϵ2
- L1-norm, v ⟵ v / ∥ v ∥ 1 + ϵ \mathrm{v \longleftarrow v / {\|v\|_{1} + \epsilon}} v⟵v/∥v∥1+ϵ
- L1-sqrt, v ⟵ v / ( ∥ v ∥ 1 + ϵ ) \mathrm{v \longleftarrow \sqrt{v/(\|v\|_{1} + \epsilon})} v⟵v/(∥v∥1+ϵ)
- L2-Hys, 它可以通过先进行L2-norm，对结果进行截短（clipping）（即v的最大值被限制为0.2之内），然后再重新归一化得到。
采用L2-Hys，L2-norm 和 L1-sqrt方式所取得的效果是一样的，L1-norm稍微表现出一点点不可靠性。但是对于没有被归一化的数据来说，这四种方法都表现出来显着的改进。参考一下matlab官方提供的代码
% -------------------------------------------------------------------------
% Normalize vector using L2-Hys
% -------------------------------------------------------------------------
function x = normalizeL2Hys(x)
classToUse = class(x);
x = x./(norm(x,2) + eps(classToUse)); % L2 norm, NORM(X,2) returns the 2-norm of X.
x(x > 0.2) = 0.2; % Clip to 0.2
x = x./(norm(x,2) + eps(classToUse)); % repeat L2 norm
(6-4) 区间(块)有两个主要的几何形状——矩形区间（R-HOG）和环形区间（C-HOG）(跟踪中没用，省略)。
R-HOG区间（blocks）：大体上是一些方形的格子，它可以有三个参数来表征：每个区间中细胞单元的数目、每个细胞单元中像素点的数目、每个细胞的直方图通道数目。例如：行人检测的最佳参数设置是：3×3细胞/区间、6×6像素/细胞、9个直方图通道。则一块的特征数为：3*3*9；作者还发现，对于R-HOG，在对直方图做处理之前，给每个区间（block）加一个高斯空域窗口（Gaussian spatial window）是非常必要的，因为这样可以降低边缘的周围像素点（pixels around the edge）的权重。R-HOG是各区间被组合起来用于对空域信息进行编码（are used in conjunction to encode spatial form information）。
(6-5) HOG描述符（不同于OpenCV定义）：我们将归一化之后的块描述符（向量）就称之为HOG描述符。
(6-6) 块划分带来的问题：块与块之间是相互独立的吗？
答：通常的将某个变量范围固定划分为几个区域，由于边界变量与相邻区域也有相关性，所以变量只对一个区域进行投影而对相邻区域完全无关时会对其他区域产生混叠效应。
(6-7) 生成HOG特征描述向量
将所有“block”的HOG描述符组合在一起，形成最终的feature vector，该feature vector就描述了detect window的图像内容。
3 问题回答
- 为什么使用方向梯度直方图？
Capture local shape information - 为什么使用 “cell”？
achieve a small amount of spatial invariance - 为什么使用“overlapping blocks”？
在众多的local contrast normalisation方法中，采用overlapping blocks得到的效果最好。
4 代码
4.1 matlab 代码
[注]：代码并非按论文来，只是一个简化版本。可以参考matlab自带的代码，更容易理解和学习。
function [hog_Feature] = HOGExtract(img, options)
% HOGEXTRACT extract Histogram of oriented gradient
% ------------------------------------------------------------------------------
% img - Input image.
% options - struct of para.
% cellH, cellW: cell Size.
% blockH, blockW: block Size.
% winH, winW: boundingbox Size.
% stride: block sliding step Size.
% bins: length of histogram
% flag: Gradient direction mapping interval
% (0: [0, pi], 1: [0, 2 * pi])
% epsilon: Normalization factor
% @hog_Feature - Result feature vector, Size: 1 * M
% M = ( (winH - blockH) / stride + 1 ) * ( (winW - blockW) / stride + 1 ) ...
% * (blockW / cellW) * (blockH / cellH) * bins
%
% Extract Step:
% Step 1: Rgb to Gray
% Step 1.1: Normalise Gamma (Optional)
% Step 2: Compute gadient, operator: [-1 0 1] [-1 0 1]'
%
% ||grad|| = |grad_x| + |grad_y| or sqrt(grad_x^2 + grad_y^2)
% gradOri = atan(grad_y / grad_x), Range: [-pi / 2, pi/2]
%
% If flag == 0, gradOri -> [0 pi]
% If flag == 1, gradOri -> [-pi/2 pi/2]
% |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| %
% gradOri(i, j) = gradOri(i,j) < 0 ? gradOri(i, j) + pi, gradOri(i,j)
% if flag == 1 gradOri -> [0 2*pi]
% 1. grad_x >= 0 && grad_y >= 0 gradOri(i,j) = atan(grad_y / grad_x)
% 2. grad_x < 0 && grad_y >= 0 gradOri(i,j) = atan(grad_y / grad_x) + pi
% 3. grad_x < 0 && grad_y < 0 gradOri(i,j) = atan(grad_y / grad_x) + pi
% 4. grad_x >= 0 && grad_y < 0 gradOri(i,j) = atan(grad_y / grad_x) + 2*pi
% |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| %
% Step 3: 4 layer loop
% layer 1 and layer 2: blockW and blockH
% layer 3 and layer 4: cellW and cellW
%
% 1) Compute cell historgrim: Using trilinear interpolation compute.
% 2) Use gaussian filter decease overlap.
% 3) Use L2-norm.
%
% Step 4: Combine the eigenvectors of all blocks
%
% Author: L.L.He
% Revised: Animate Deng
% Date: 2018-11-21
%
% ------------------------------------------------------------------------------
tic;
assert(nargin >= 1, 'Please input a img matrix!');
if ~exist('options', 'var')
% Default value
options = struct;
options.cellH = 8;
options.cellW = 8;
options.blockH = 16;
options.blockW = 16;
options.winH = 64;
options.winW = 128;
options.stride = 8;
options.bins = 9;
options.flag = 1;
options.epsilon = 1e-4;
end
% Process the input image
[r, c, d] = size(img);
if d ~= 1 % If img is rgb, rgb --> gray
img = rgb2gray(img);
end
img = double(img); % Ensure accuracy
if r ~= options.winH && c ~= options.winW
% Resize img by bilinear interpolate, --> [128 64]
img = imresize(img, [options.winH options.winW], 'bilinear');
end
% Step2: Use opeator [-1 0 1] and -[-1 0 1]' compute gradient
% ---------------------------------------------------------------------------- %
mask = [-1 0 1];
[grad, gradOri] = calGrad(img, mask, options.flag);
% Compute gaussian kernel size
sigma = min(options.blockH, options.blockW) * 0.5;
Sigma = sigma .^ 2;
[X, Y] = meshgrid(0 : options.blockW-1, 0 : options.blockH-1);
X = X - (options.blockW - 1) / 2;
Y = Y - (options.blockH - 1) / 2;
gaussWeight = exp(-(X.^2 + Y.^2) / (2 * Sigma));
% Compute Historgrim
r_size = (options.winH - options.blockH) / options.stride + 1;
c_size = (options.winW - options.blockW) / options.stride + 1;
b_size = options.bins * (options.blockH * options.blockW) / ...
(options.cellH * options.cellW);
blockHist = zeros(r_size, c_size, b_size);
% Step 3: Compute Hist
% ---------------------------------------------------------------------------- %
for i = 1 : options.stride : (options.winH - options.blockH + 1)
for j = 1 : options.stride : (options.winW - options.blockW + 1)
block_grad = grad(i : i+options.blockH-1, j : j+options.blockW-1);
block_gradOri = gradOri(i : i+options.blockH-1, j : j+options.blockW-1);
block_r = floor(i / options.stride) + 1;
block_c = floor(j / options.stride) + 1;
blockHist(block_r, block_c, :) = calHist(block_grad .* gaussWeight, ...
block_gradOri, options);
end
end
% Step 4: Combine the eigenvectors of all blocks
% ---------------------------------------------------------------------------- %
hog_Feature = reshape(blockHist, [1 numel(blockHist)]);
toc;
end
% ================================ functions ================================= %
function [grad, gradOri] = calGrad(img, mask, flag)
% Compute gradient and angel
assert(nargin == 3);
maskX = mask;
maskY = -mask'; % Note
grad = zeros( size(img) );
gradOri = zeros( size(img) );
grad_x = imfilter(img, maskX);
grad_y = imfilter(img, maskY);
grad = sqrt(grad_x .^ 2 + grad_y .^ 2); % Compute Amplitude
if flag == 0
% Ori --> [0 pi]
gradOri = atan(grad_y ./ (grad_x + eps));
idx = find(gradOri < 0);
gradOri(idx) = gradOri(idx) + pi;
else
% Since the atan return value is in the range [-pi/2, pi/2],
% it needs to be restored to [0, 2 * pi]
idx_1 = find(grad_x >= 0 & grad_y >= 0);
gradOri(idx_1) = atan(grad_y(idx_1) ./ (grad_x(idx_1) + eps));
idx_23 = find(grad_x < 0);
gradOri(idx_23) = atan(grad_y(idx_23) ./ (grad_x(idx_23) + eps)) + pi;
idx_4 = find(grad_x >= 0 & grad_y < 0);
gradOri(idx_4) = atan(grad_y(idx_4) ./ (grad_x(idx_4) + eps)) + 2 * pi;
end
end
% ================================ functions ================================= %
function hist = calHist(block_grad, block_gradOri, options)
% Compute Hist
bins = options.bins;
cellH = options.cellH;
cellW = options.cellW;
blockH = options.blockH;
blockW = options.blockW;
assert(mod(blockH, cellH) == 0 && mod(blockW, cellW) == 0);
hist = zeros(blockH / cellH, blockW / cellW, bins); % 3 D --> [2 2 9]
% bins [0 : 20 : pi], 9 region
if options.flag == 0
anglePerBin = pi / bins;
correctVar = pi; % Correct curri < 0
else
anglePerBin = 2 * pi / bins; % 40
correctVar = 2 * pi; % 360
end
halfAngle = anglePerBin / 2; % 20
for i = 1 : blockH
for j = 1 : blockW
cell_r = floor( (i-1) / cellH ) + 1;
cell_c = floor( (j-1) / cellW ) + 1;
% Calculate the two bins connected to the current pixel and vote
% ======================== E.g ========================== %
% if : block_gradOri(i, j) = 165
% So : The range in which voting is considered is [160 0] [0 20]
% Compute weight and hist
% ======================================================= %
currOri = block_gradOri(i, j) - halfAngle;
% In order to connect the first bin to the last bin,
% it is equivalent to 0 degrees and 180(or 360) degrees.
if currOri <= 0
currOri = currOri + correctVar;
end
% Compute the point loc.
pre_idxOfbin = floor(currOri / anglePerBin) + 1;
pro_idxOfbin = mod(pre_idxOfbin, bins) + 1;
% Vote for two adjacent bins (distance to center point as weight)
center = (2 * pre_idxOfbin-1) * halfAngle;
dist_w = (currOri + halfAngle - center) / anglePerBin;
hist(cell_r, cell_c, pre_idxOfbin) = hist(cell_r, cell_c, ...
pre_idxOfbin) + (1 - dist_w) * block_grad(i, j);
hist(cell_r, cell_c, pro_idxOfbin) = hist(cell_r, cell_c, ...
pro_idxOfbin) + dist_w * block_grad(i, j);
end
end
hist = reshape(hist, [1 numel(hist)]);
% L2-norm
hist = hist ./ sqrt(hist * hist' + options.epsilon .^ 2);
end
4.2 OpenCV 代码
/* 待补充 */
5 测试
clc; clear; close all;
img = imread('C:\Users\xx\Pictures\Saved Pictures\matlab\house.png');
[hog_Feature0] = HOGExtract(img);
tic;
[hog_Feature1, validPoints] = extractHOGFeatures(img);
toc;
figure, imshow(img, []);
hold on;
plot(validPoints);
[M, N] = size(hog_Feature0);
扫一扫
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
