OpenCV3特征提取与目标检测之HOG（一）——HOG的概述与原理

1.HOG(Histogram of Oriented Gradient）是方向梯度直方图的意思，是一种特性描述子。通过计算与统计图像局部区域的梯度方向直方图来构成特征。边缘是图像颜色剧变的区域，在一副图像中，局部目标的表象与形状能够被梯度或边缘的方向密度分布很好地描述，而梯度主要存在于存在于局部目标边缘的地方。
(1) 局部目标的边缘，可以把图转为灰度图后按二值映射到0和1输出，就很明显的看出来。

而梯度可分解为 x 方向的梯度 G{x} 和 y 方向的梯度 G{y} 。某个像素点的 x 方向的梯度的计算可以通过这个像素点左右两边的像素值的差值的绝对值计算出来，而 y 方向的梯度可以通过该像素点上下两边的像素值的差值的绝对值计算。而根据下面的两个公式可以计算每一个像素点的梯度方向和梯度幅值。

(2)直方图是图像中像素强度分布的图形表达方式，一张图像可以看成一个数字矩阵，矩阵的每个元素取值在0~255之间，而用直方图来统计了每个强度值所具有的像素个数。那么已知数字的范围包含256个值，将这个范围分割成子区域就是bins，然后再统计在每一个bin的像素数目。这样统计右边的数字矩阵我们可以得到左图(x轴表示bin, y轴表示各个bin中的像素个数)。

（3） HOG是通过上面公式计算出来的梯度方向的角度是一个范围在0-360度的弧度值，为了计算简单，将梯度向的范围约束为0-180度，并且分割为9个方向，每个方向20度，再将约束后的角度除以20，则现在的梯度方向角度值就变为范围在[0,9)。

2.HOG特征提取具体实现方法是将图像分成小的连通区域，叫细胞单元（cell），将每个小Cell里面的梯度幅值按照9个方向进行统计，计算完之后，将会产生一个横坐标X为梯度方向，纵坐标Y为梯度幅值的方向梯度直方图。然后采集细胞单元中各像素点的梯度的或边缘的方向直方图,最后把这些直方图组合起来就可以构成特征描述器。把这些局部直方图在图像的更大的范围内（block）进行对比度归一化。 归一化是为了克服光照不均匀的变化以及前景和背景测对比差异。
大概步骤：
（1）灰度化。
（2）对输入图像进行颜色空间的标准化（归一化）；调节图像的对比度，降低图像局部的阴影和光照变化所造成的影响，同时可以抑制噪音的干扰。
（3）计算图像每个像素的梯度（包括大小和方向）；主要是为了捕获轮廓信息，同时进一步弱化光照的干扰。
（4）将图像划分成小cells（例如8*8像素/cell）。
（5）统计每个cell的梯度直方图（不同梯度的个数），即可形成每个cell的特性描述子。
（6）将每几个cell组成一个block（2X2个cell/block），一个block内所有cell的特性描述子串联起来便得到该block的HOG特性描述子。
（7）将图像内的所有block的HOG特征描述子串联起来就可以得到该图像的HOG特征描述子。这个就是最终的可供分类使用的特征向量了。

3.cell、block、windowsSize、stride的关系。

上图中单个cell的为8X8个像素，把cell对应的方向直方图转换为单维向量，按规定组距对对应方向梯度个数进行编码，得到单个cell的9个特征，每个block包含2X2个cell,那么每个block包含2X2个cell也就是2X2X9=36个特征，而每个block移动（stride）这里选择overlap，就是为2分之一重叠，一个64X128大小的图像横着有15个block,坚着有7个，最后得到的特征数为36X7X15=3780维。
5.代码实现
这时的代码环境是win7 64位，vs2015,opencv3.3。
（1）代码流程

（2）具体代码
头文件

#pragma once
#include <opencv2\core\core.hpp>
#include <opencv2\highgui.hpp>
#include <opencv2\xfeatures2d\nonfree.hpp>
#include <opencv2\opencv.hpp>
#include <iostream>

using namespace std;
using namespace cv;

class HogFeatureExtraction
{
public:
	//构造函数
	HogFeatureExtraction();
	//构造函数
	HogFeatureExtraction(int _bins, int _theta, int _cell_size, int _block_size, float _R);
	~HogFeatureExtraction();

	//计算积分图
	vector<Mat> calculateIntegralHOG(Mat& srcMat);
	//单个cell的HOG
	void cacHOGinCell(Mat& m_HOGCell, Rect roi, vector<Mat>& m_integrals);
	//单个block的HOG
	Mat  getBlockHog(cv::Point pt, std::vector<cv::Mat>& integrals);
	//整个图像的HOG
	Mat cacHOGFeature(cv::Mat srcImage);

private:
	int bins;
	int theta;
	int cell_size;
	int block_size;
	float R;
};

实现文件

#include "HogFeatureExtraction.h"

HogFeatureExtraction::HogFeatureExtraction()
{
	bins = 9;
	theta = 180 / bins;
	cell_size = 20;
	block_size = 2;
	R = cell_size*(block_size)*0.5;
}
HogFeatureExtraction::HogFeatureExtraction(int _bins, int _theta, int _cell_size, int _block_size)
{
	bins = _bins;
	theta = 180 / bins;
	block_size = _block_size;
	cell_size = _cell_size;
	R = _cell_size*(_block_size)*0.5;
}

HogFeatureExtraction::~HogFeatureExtraction()
{
}
// 计算积分图
vector<Mat> HogFeatureExtraction::calculateIntegralHOG(Mat& m_src)
{
	Mat m_sobel_x,m_sobel_y,m_magn,m_angle;

	//x 方向上的差分阶数
	Sobel(m_src, m_sobel_x, CV_32F, 1, 0);
	//y 方向上的差分阶数
	Sobel(m_src, m_sobel_y, CV_32F, 0, 1);

	//根据每一个点X方向和Y方向上的梯度，实现笛卡尔坐标和极坐标的转换
	//得到梯度幅值和梯度方向
	cartToPolar(m_sobel_x, m_sobel_y, m_magn, m_angle, true);
	//将笛卡尔坐标转换为极坐标之后，角度的范围在[0,360],但要转成[0,180]
	//如果m_angle<0，则加180
	add(m_angle, Scalar(180), m_angle, m_angle<0);
	//如果m_angle>=180，则减180
	add(m_angle, Scalar(-180), m_angle, m_angle >= 180);

	//将角度矩阵转换为一个灰度值范围在0~9之间的图像
	m_angle = m_angle /theta;
	//新建9个矩阵
	vector<Mat> m_bins(bins);
	for (int i = 0; i < bins; i++)
	{
		//初始化为全0矩阵
		m_bins[i] = Mat::zeros(m_src.size(), CV_32F);
	}	
	//把图像的梯度幅值矩阵按九个不同方向的梯度角度，将每个角度范围内相应点的梯度幅值存储在相应的矩阵图像之上
	for (int y = 0; y < m_src.rows; y++)
	{
		for (int x = 0; x < m_src.cols; x++)
		{
			int ind = m_angle.at<float>(y, x);
			m_bins[ind].at<float>(y, x) += m_magn.at<float>(y, x);
		}
	}
	/*根据上面生成的9张不同角度的梯度幅值矩阵生成9张不同的梯度幅值的积分图像，至此以后，
	  积分图像的每一点就代表，这一点左上角，所有梯度幅值之和；生成的9幅积分图也就是9个
	  bins，不同bins上的HOG强度*/
	vector<Mat> integrals(bins);
	for (int i = 0; i < bins; i++)
	{
		//得到积分图
		integral(m_bins[i], integrals[i]);
	}
	return integrals;
}
//计算单个cell的HOG特征
void HogFeatureExtraction::cacHOGinCell(Mat& m_HOGCell, Rect roi, vector<Mat>& vm_integrals)
{
	//通过9幅积分图像快速实现HOG的计算，HOG这个直方图有9个bins，每个bins就对应一张积分图像
	//确定单个矩形cell的左上角点坐标
	int x0 = roi.x;                             
	int y0 = roi.y;
	//确定单个矩形cell的右下角点坐标
	int x1 = x0 + roi.width;
	int y1 = y0 + roi.height;                   

	for (int i = 0; i <bins; i++)
	{
		//初始化容器
		Mat integral = vm_integrals[i];
		float a = integral.at<double>(y0, x0);
		float b = integral.at<double>(y1, x1);
		float c = integral.at<double>(y0, x1);
		float d = integral.at<double>(y1, x0);
		/*每循环一次，计算一个梯度方向上的HOG特征，*/
		/*每循环一次，就计算梯度方向直方图上的一个bins*/
		m_HOGCell.at<float>(0, i) = b - c - d + a;
	}
}
//计算单个Block的HOG梯度方向直方图
Mat HogFeatureExtraction::getBlockHog(Point pt, vector<Mat>& m_integrals)
{
	if (pt.x - R<0 || pt.y - R<0 || pt.x + R >= m_integrals[0].cols || pt.y + R >= m_integrals[0].rows)
	{
		return Mat();
	}
	//初始化一个矩阵
	Mat hist(Size(bins*block_size*block_size, 1), CV_32F);
	Point  t1(0, pt.y - R);
	int c = 0;
	//遍历单个block,就是遍历4个cell，并且将4个cell的HOG特征向量组成了一个维数比较大的block的HOG特征向量
	for (int i = 0; i<block_size; i++)
	{
		t1.x = pt.x - R;
		for (int j = 0; j<block_size; j++)
		{
			//获取当前窗口，进行局部HOG直方图计算
			Rect roi(t1, t1 + Point(cell_size, cell_size));
			Mat  hist_temp = hist.colRange(c, c + bins);
			//根据roi确定的矩形区域，计算单个cell的HOG直方图(其本质就是一个行特征向量)
			cacHOGinCell(hist_temp, roi, m_integrals);
			t1.x += cell_size;
			c += bins;
		}
		t1.y = cell_size;
	}
	 //归一化
	normalize(hist, hist, 1, 0, NORM_L2);
	return hist;
}

//计算整幅图像的HOG梯度方向直方图
Mat HogFeatureExtraction::cacHOGFeature(Mat src_image)
{
	Mat gray_image, dst_image;
	vector<Mat> wm_HOG;
	//灰度图
	cvtColor(src_image, gray_image, CV_RGB2GRAY);
	//转换类型
	gray_image.convertTo(gray_image, CV_8UC1);
	//生成9个不同梯度方向上的梯度幅值的积分图像
	vector<Mat> integrals = calculateIntegralHOG(gray_image);
	dst_image = gray_image.clone();
	dst_image = dst_image * 0.5;
	
	Mat m_HOGBlock(Size(bins, 1), CV_32F);
	//遍历全图像，计算最终的梯度方向直方图HOG
	for (int y = cell_size / 2; y < gray_image.rows; y += cell_size)
	{
		for (int x = cell_size / 2; x < gray_image.cols; x += cell_size)
		{
			//获取当前block的HOG特征，每个block由四个clee组成，每个cel由8*8个像素组成
			Mat hist = getBlockHog(Point(x, y), integrals);
			if (hist.empty())
			{
				continue;
			}
			//赋值为全0的矩阵
			m_HOGBlock = Scalar(0);
			//计算的就是单个Block的梯度方向直方图HOG
			for (int i = 0; i < bins; i++)
			{
				for (int j = 0; j < block_size; j++)
				{
					m_HOGBlock.at<float>(0, i) += hist.at<float>(0, i + j*bins);
				}
			}
			//对其得到的每个Block的的矩阵进行归一化，使其转变为一个block的HOG特征向量
			normalize(m_HOGBlock, m_HOGBlock, 1, 0, CV_L2);
		
			//每得到一个Block的HOG特征向量就存入wm_HOG，等到整个图像的HOG特征向量
			wm_HOG.push_back(m_HOGBlock);
			Point center(x, y);
			//绘制HOG特征图
			for (int i = 0; i < bins; i++)
			{
				double d_theta = (i * theta) * CV_PI / 180.0;
				Point rd(cell_size*0.5*cos(d_theta), cell_size*0.5*sin(d_theta));
				Point rp = center - rd;
				Point lp = center + rd;
				line(dst_image, rp, lp, Scalar(255 * m_HOGBlock.at<float>(0, i), 255, 255));
			}
		}
	}
	return dst_image;
}

主函数

#include "HogFeatureExtraction.h"

using namespace std;
using namespace cv;
using namespace cv::ml;

int main(void)
{
	Mat src = imread("C:/Users/matt/Desktop/demo/04.jpg");
	namedWindow("src", WINDOW_NORMAL);
	imshow("src", src);
	HogFeatureExtraction hogFeature(9, 180, 8, 2);
	Mat dst = hogFeature.cacHOGFeature(src);
	namedWindow("HOG", WINDOW_NORMAL);
	imshow("HOG", dst);
	waitKey(0);
	return 0;
}

运行效果

这个像素只有128*64，效果不明显，可以放张大一些的图像试试。
结语：
1.以上只是介绍了HOG特性的提取，之后会试着跑官方行人检测的demo,和如何应用HOG训练自己的分类模型。
2.关于整个工程的源码，运行程序时的bug,或者有如何优代码的想法都可以加这个群(487350510)互相讨论学习。