本文原文链接：https://blog.csdn.net/ppp8300885/article/details/71078555
本文用80行代码的Python实现了HOG算法，代码在Github Hog-feature，虽然OpenCV有实现好的Hog描述器算法，但是本文目的是完全理解HOG特征提取的具体方法和实现原理，以及检验相关参数对实验结果的影响，提升检测到的特征的性能以及优化代码的运行速度。
1. 方法简介
方向梯度直方图（Histogram of Oriented Gradient, HOG）特征是一种在计算机视觉和图像处理中用来进行物体检测的描述子。通过计算和统计局部区域的梯度方向直方图来构成特征。Hog特征结合SVM分类器已经被广泛应用于图像识别中，尤其在行人检测中获得了极大的成功。现如今如今虽然有很多行人检测算法不断提出，但基本都是以HOG+SVM的思路为主。
主要思想：在一幅图像中，局部目标的表象和形状（appearance and shape）能够被梯度或边缘的方向密度分布很好地描述。其本质是梯度的统计信息，而梯度主要存在于边缘所在的地方。
实现过程：简单来说，首先需要将图像分成小的连通区域，称之为细胞单元。然后采集细胞单元中各像素点的梯度或边缘的方向直方图。最后把这些直方图组合起来就可以构成特征描述器。
算法优点：与其他的特征描述方法相比，HOG有较多优点。由于HOG是在图像的局部方格单元上操作，所以它对图像几何的和光学的形变都能保持很好的不变性，这两种形变只会出现在更大的空间领域上。其次，在粗的空域抽样、精细的方向抽样以及较强的局部光学归一化等条件下，只要行人大体上能够保持直立的姿势，可以容许行人有一些细微的肢体动作，这些细微的动作可以被忽略而不影响检测效果。因此HOG特征是特别适合于做图像中的人体检测的。
2. HOG流程
HOG特征提取算法的整个实现过程大致如下：
读入所需要的检测目标即输入的image
将图像进行灰度化（将输入的彩色的图像的r,g,b值通过特定公式转换为灰度值）
采用Gamma校正法对输入图像进行颜色空间的标准化（归一化）
计算图像每个像素的梯度（包括大小和方向），捕获轮廓信息
统计每个cell的梯度直方图（不同梯度的个数），形成每个cell的descriptor
将每几个cell组成一个block（以3*3为例），一个block内所有cell的特征串联起来得到该block的HOG特征descriptor
将图像image内所有block的HOG特征descriptor串联起来得到该image（检测目标）的HOG特征descriptor，这就是最终分类的特征向量
HOG参数设置是：2*2细胞／区间、8*8像素／细胞、8个直方图通道,步长为1。
特征提取流程图如下:
3. Python实现
3.1 数据准备
读入彩色图像，并转换为灰度值图像, 获得图像的宽和高。采用Gamma校正法对输入图像进行颜色空间的标准化（归一化），目的是调节图像的对比度，降低图像局部的阴影和光照变化所造成的影响，同时可以抑制噪音。采用的gamma值为0.5。
#first part
#first part import cv2 import numpy as np img = cv2.imread('person_037.png', cv2.IMREAD_GRAYSCALE) # cv2.imshow('Image', img) # cv2.imwrite("Image-test.jpg", img) # cv2.waitKey(0) img = np.sqrt(img / float(np.max(img))) # cv2.imshow('Image', img) # cv2.imwrite("Image-test2.jpg", img) # cv2.waitKey(0)
3.2 计算每个像素的梯度
计算图像横坐标和纵坐标方向的梯度，并据此计算每个像素位置的梯度方向值；求导操作不仅能够捕获轮廓，人影和一些纹理信息，还能进一步弱化光照的影响。在求出输入图像中像素点（x,y）处的水平方向梯度、垂直方向梯度和像素值，从而求出梯度幅值和方向。
常用的方法是：首先用[-1,0,1]梯度算子对原图像做卷积运算，得到x方向（水平方向，以向右为正方向）的梯度分量gradscalx，然后用[1,0,-1]T梯度算子对原图像做卷积运算，得到y方向（竖直方向，以向上为正方向）的梯度分量gradscaly。然后再用以上公式计算该像素点的梯度大小和方向。
# second part
# second part height, width = img.shape gradient_values_x = cv2.Sobel(img, cv2.CV_64F, 1, 0, ksize=5) gradient_values_y = cv2.Sobel(img, cv2.CV_64F, 0, 1, ksize=5) gradient_magnitude = cv2.addWeighted(gradient_values_x, 0.5, gradient_values_y, 0.5, 0) gradient_angle = cv2.phase(gradient_values_x, gradient_values_y, angleInDegrees=True) print gradient_magnitude.shape, gradient_angle.shape
Out
(640, 480) (640, 480)
1
3.3 为每个细胞单元构建梯度方向直方图
我们将图像分成若干个“单元格cell”，默认我们将cell设为8*8个像素。假设我们采用8个bin的直方图来统计这6*6个像素的梯度信息。也就是将cell的梯度方向360度分成8个方向块，例如：如果这个像素的梯度方向是0-22.5度，直方图第1个bin的计数就加一，这样，对cell内每个像素用梯度方向在直方图中进行加权投影（映射到固定的角度范围），就可以得到这个cell的梯度方向直方图了，就是该cell对应的8维特征向量而梯度大小作为投影的权值。
# third part
# third part cell_size = 8 bin_size = 8 angle_unit = 360 / bin_size gradient_magnitude = abs(gradient_magnitude) cell_gradient_vector = np.zeros((height / cell_size, width / cell_size, bin_size)) print cell_gradient_vector.shape def cell_gradient(cell_magnitude, cell_angle): orientation_centers = [0] * bin_size for k in range(cell_magnitude.shape[0]): for l in range(cell_magnitude.shape[1]): gradient_strength = cell_magnitude[k][l] gradient_angle = cell_angle[k][l] min_angle = int(gradient_angle / angle_unit)%8 max_angle = (min_angle + 1) % bin_size mod = gradient_angle % angle_unit orientation_centers[min_angle] += (gradient_strength * (1 - (mod / angle_unit))) orientation_centers[max_angle] += (gradient_strength * (mod / angle_unit)) return orientation_centers for i in range(cell_gradient_vector.shape[0]): for j in range(cell_gradient_vector.shape[1]): cell_magnitude = gradient_magnitude[i * cell_size:(i + 1) * cell_size, j * cell_size:(j + 1) * cell_size] cell_angle = gradient_angle[i * cell_size:(i + 1) * cell_size, j * cell_size:(j + 1) * cell_size] print cell_angle.max() cell_gradient_vector[i][j] = cell_gradient(cell_magnitude, cell_angle)
3.4 可视化Cell梯度直方图
将得到的每个cell的梯度方向直方图绘出，得到特征图
# fourth part
# fourth part import math import matplotlib.pyplot as plt hog_image= np.zeros([height, width]) cell_gradient = cell_gradient_vector cell_width = cell_size / 2 max_mag = np.array(cell_gradient).max() for x in range(cell_gradient.shape[0]): for y in range(cell_gradient.shape[1]): cell_grad = cell_gradient[x][y] cell_grad /= max_mag angle = 0 angle_gap = angle_unit for magnitude in cell_grad: angle_radian = math.radians(angle) x1 = int(x * cell_size + magnitude * cell_width * math.cos(angle_radian)) y1 = int(y * cell_size + magnitude * cell_width * math.sin(angle_radian)) x2 = int(x * cell_size - magnitude * cell_width * math.cos(angle_radian)) y2 = int(y * cell_size - magnitude * cell_width * math.sin(angle_radian)) cv2.line(hog_image, (y1, x1), (y2, x2), int(255 * math.sqrt(magnitude))) angle += angle_gap plt.imshow(hog_image, cmap=plt.cm.gray) plt.show()
out
3.5 统计Block的梯度信息
把细胞单元组合成大的块(block），块内归一化梯度直方图
由于局部光照的变化以及前景-背景对比度的变化，使得梯度强度的变化范围非常大。这就需要对梯度强度做归一化。归一化能够进一步地对光照、阴影和边缘进行压缩。
把各个细胞单元组合成大的、空间上连通的区间（blocks）。这样，一个block内所有cell的特征向量串联起来便得到该block的HOG特征。这些区间是互有重叠的，
本次实验采用的是矩阵形区间，它可以有三个参数来表征：每个区间中细胞单元的数目、每个细胞单元中像素点的数目、每个细胞的直方图通道数目。
本次实验中我们采用的参数设置是：2*2细胞／区间、8*8像素／细胞、8个直方图通道,步长为1。则一块的特征数为2*2*8。
# fifth part
# fifth part hog_vector = [] for i in range(cell_gradient_vector.shape[0] - 1): for j in range(cell_gradient_vector.shape[1] - 1): block_vector = [] block_vector.extend(cell_gradient_vector[i][j]) block_vector.extend(cell_gradient_vector[i][j + 1]) block_vector.extend(cell_gradient_vector[i + 1][j]) block_vector.extend(cell_gradient_vector[i + 1][j + 1]) mag = lambda vector: math.sqrt(sum(i ** 2 for i in vector)) magnitude = mag(block_vector) if magnitude != 0: normalize = lambda block_vector, magnitude: [element / magnitude for element in block_vector] block_vector = normalize(block_vector, magnitude) hog_vector.append(block_vector) print np.array(hog_vector).shape
out
(4661, 32)
1
共有4661个block，每个block有32维的特征
4. 80行代码封装
通过建类和封装函数，最终代码如下：
import cv2 import numpy as np import math import matplotlib.pyplot as plt class Hog_descriptor(): def __init__(self, img, cell_size=16, bin_size=8): self.img = img self.img = np.sqrt(img / np.max(img)) self.img = img * 255 self.cell_size = cell_size self.bin_size = bin_size self.angle_unit = 360 / self.bin_size assert type(self.bin_size) == int, "bin_size should be integer," assert type(self.cell_size) == int, "cell_size should be integer," assert type(self.angle_unit) == int, "bin_size should be divisible by 360" def extract(self): height, width = self.img.shape gradient_magnitude, gradient_angle = self.global_gradient() gradient_magnitude = abs(gradient_magnitude) cell_gradient_vector = np.zeros((height / self.cell_size, width / self.cell_size, self.bin_size)) for i in range(cell_gradient_vector.shape[0]): for j in range(cell_gradient_vector.shape[1]): cell_magnitude = gradient_magnitude[i * self.cell_size:(i + 1) * self.cell_size, j * self.cell_size:(j + 1) * self.cell_size] cell_angle = gradient_angle[i * self.cell_size:(i + 1) * self.cell_size, j * self.cell_size:(j + 1) * self.cell_size] cell_gradient_vector[i][j] = self.cell_gradient(cell_magnitude, cell_angle) hog_image = self.render_gradient(np.zeros([height, width]), cell_gradient_vector) hog_vector = [] for i in range(cell_gradient_vector.shape[0] - 1): for j in range(cell_gradient_vector.shape[1] - 1): block_vector = [] block_vector.extend(cell_gradient_vector[i][j]) block_vector.extend(cell_gradient_vector[i][j + 1]) block_vector.extend(cell_gradient_vector[i + 1][j]) block_vector.extend(cell_gradient_vector[i + 1][j + 1]) mag = lambda vector: math.sqrt(sum(i ** 2 for i in vector)) magnitude = mag(block_vector) if magnitude != 0: normalize = lambda block_vector, magnitude: [element / magnitude for element in block_vector] block_vector = normalize(block_vector, magnitude) hog_vector.append(block_vector) return hog_vector, hog_image def global_gradient(self): gradient_values_x = cv2.Sobel(self.img, cv2.CV_64F, 1, 0, ksize=5) gradient_values_y = cv2.Sobel(self.img, cv2.CV_64F, 0, 1, ksize=5) gradient_magnitude = cv2.addWeighted(gradient_values_x, 0.5, gradient_values_y, 0.5, 0) gradient_angle = cv2.phase(gradient_values_x, gradient_values_y, angleInDegrees=True) return gradient_magnitude, gradient_angle def cell_gradient(self, cell_magnitude, cell_angle): orientation_centers = [0] * self.bin_size for i in range(cell_magnitude.shape[0]): for j in range(cell_magnitude.shape[1]): gradient_strength = cell_magnitude[i][j] gradient_angle = cell_angle[i][j] min_angle, max_angle, mod = self.get_closest_bins(gradient_angle) orientation_centers[min_angle] += (gradient_strength * (1 - (mod / self.angle_unit))) orientation_centers[max_angle] += (gradient_strength * (mod / self.angle_unit)) return orientation_centers def get_closest_bins(self, gradient_angle): idx = int(gradient_angle / self.angle_unit) mod = gradient_angle % self.angle_unit return idx, (idx + 1) % self.bin_size, mod def render_gradient(self, image, cell_gradient): cell_width = self.cell_size / 2 max_mag = np.array(cell_gradient).max() for x in range(cell_gradient.shape[0]): for y in range(cell_gradient.shape[1]): cell_grad = cell_gradient[x][y] cell_grad /= max_mag angle = 0 angle_gap = self.angle_unit for magnitude in cell_grad: angle_radian = math.radians(angle) x1 = int(x * self.cell_size + magnitude * cell_width * math.cos(angle_radian)) y1 = int(y * self.cell_size + magnitude * cell_width * math.sin(angle_radian)) x2 = int(x * self.cell_size - magnitude * cell_width * math.cos(angle_radian)) y2 = int(y * self.cell_size - magnitude * cell_width * math.sin(angle_radian)) cv2.line(image, (y1, x1), (y2, x2), int(255 * math.sqrt(magnitude))) angle += angle_gap return image img = cv2.imread('person_037.png', cv2.IMREAD_GRAYSCALE) hog = Hog_descriptor(img, cell_size=8, bin_size=8) vector, image = hog.extract() print np.array(vector).shape plt.imshow(image, cmap=plt.cm.gray) plt.show()
5. 结果分析
本文最终单幅图像HOG特征的求取平均时间为1.8秒，相比最初版本所需的5.4秒有个长足的改进。
相比初期的版本hog梯度特征图
可见最终版本中
能够更加有效的区分梯度显示边缘。这是因为对各个像素的梯度进行了全局归一化，并且在描绘梯度方向时加入了梯度量级的非线性映射，使得梯度方向产生明显的深浅和长度差异，更易于区分边缘，凸显明显的梯度变化。
此外在输入图像时，采用Gamma校正对输入图像进行颜色空间的标准化能够抑制噪声，使得产生的边缘更加明显，清晰。
此外改变cell的大小和直方图方向通道的效果如下：
cell_size = 10 即 16*16个像素
可以看出增大cell的size得到的特征图更加注重基本轮廓和边缘，而忽略一些细节，某种程度上降低了噪声。
当通道数目为16个方向
梯度特征图像的细节变得更加明显，方向更多。
6. 在人脸识别，物体检测中的应用
在提取完图像的HOG特征之后，可以使用SVM进行分类训练，能完成行人检测等任务。
未来工作可参考Github的行人检测项目https://github.com/icsfy/Pedestrian_Detection