来自http://hi.baidu.com/timehandle/home
去年我做一个作业和调研的时候看了HOG的那篇论文，没怎么看明白，于是就求助于opencv的代码，结果是也没看明白。在网上做了一些求助后只好自己硬着头皮慢慢抠了两周，最后写了matlab的程序《Histograms of Oriented Gradients (HOG)特征 MATLAB 计算》，我的代码长的比较直白，但是干活不给力，很慢。大家问了一些问题，有一些比较类似，我总结了下将回答写成了这个文档。虽然我很想用图把问题描述清楚，但是可是用画图板和PHOTOSHOP三脚猫的功夫只能到这种程度了。还有就是，我不是做行人检测方向的。。。我只是碰巧碰了下HOG。。。
代码中处理globalinterpolate的情况是没有理解HOG的情况下写的，比原始HOG的想法简陋。Localinterpolate的情况是按照原始HOG实现的，是一种naïve实现方式.
关于计算梯度方向角的：
首先用[-1,0,1]梯度算子对原图像做卷积运算，得到x方向（水平方向，以向右为正方向）的梯度分量gradscalx，然后用[1,0,-1]’梯度算子对原图像做卷积运算，得到y方向（竖直方向，以向上为正方向）的梯度分量gradscaly。然后当gradscalx>=0, gradscaly>=0时，说明梯度方向是朝向第一象限的，当gradscalx>=0, gradscaly<0时，说明梯度方向是朝向第二象限的，诸如此类，结合象限信息，就可以利用反正切函数atan求出在signed和unsigned各自情况下正确的梯度角度.
关于扫描循环（四层for循环…有没有快一点的？有！但是我功力不够。。当时没编出来，就只好还是来四层for）：
假设检测窗为64(列)*128(行)大小，block为16*16大小，每个block划分为4个cell，block每次滑动8个像素(也就是一个cell的宽)，以及梯度方向划分为9个区间，在0~180度范围内统计，以下的说明都以上述假设为例.
btly与btlx分别表示block所在位置左上角点处的坐标。对于前述假设，一个检测窗内会有105个block存在，因此第一个block左上角的坐标是(1,1),第二个是(9,1)…,此行最后一个是block的左上角坐标是(49,1),然后下一个block就需要向下滑动8个像素，并回到最左边,此时的block左上角坐标为(1,9),接着block重新开始新的横向滑动…如此这般,在检测窗内最后一个block的坐标就是(49,113).
block每滑动到一个新的位置，就需要停下来计算它内部的那四个cell中的梯度方向直方图.(bj,bi)就是来存储cell左上角的坐标的（cell的坐标以block左上角为原点）.
(j,i)就表示cell中的像素在整个检测窗（64*128的图像）中的坐标.另外，我在程序里有个jorbj与iorbi，这在Localinterpolate的情况下（也就是标准的原始HOG情况），就是bj与bi.
关于hist3dbig：
这是一个三维的矩阵，用来存储三维直方图。最常见的一维的直方图是这个样子，
二维直方图呢？是这个样子，一个一个的柱子是一个统计bin，柱子的高低代表统计值的大小
三维直方图呢？是这个样子，立体的一个一个的小格子，每个小格子是一个统计bin, 小格子用来装统计值。以上面的例子，那么对一个block来说，它的直方图是下面这样的：
再来说线性插值，线性插值时，一个统计值需被“按一定比例分配”到这个统计点最邻近的区间中去，下面的图显示了一维直方图时，落在虚线标记范围内的统计点，它最近邻的区间就是标有红色圆点的两个区间
若是二维直方图，那落在如下虚线矩形中的统计点，周围的这四个统计区间就是它最近邻的区间。这个虚线矩形由四个统计区间各自的1/4组成。
三维直方图，对一个统计点来说，它的最近邻的区间有八个，如下图，可以想象一下，只有当这个统计点落在由如下八个统计区间各自的1/8组成的一个立方体内内时，这八个区间才是对统计点最近邻的。
统计时如何分配权重呢？以一维直方图简单说一下线性插值的意思，对于下面绿色小方点(x)的统计值来说，假设标红点的两个bin的中心位置分别为x1，x2，那么对于x，它的分配权重为左边bin: 1-(x-x1)/s, 即 1-a/s = b/s, 右边bin: 1-(x2-x)/s, 即1-b/s = a/s.
类似，那么对三维直方图来说，统计时的累积式（从Dalal的论文里截来的）就是：
上面，w 就是准备被分配的统计值。(x1,y1,z1)…共八个点表示八个统计区间的中心位置坐标，上式用h(x1,y1,z1)这样的标记来表示所要累积的统计区间。我在编程时就使用的这个式子，只不过我用bin的下标号来表示bin块，就像前面三维直方图示意中(binx=1,biny=2,binθ=9)，不过在程序中θ轴是用z轴表示了。
binx1 = floor((jorbj-1+cellpw/2)/cellpw) + 1;
biny1 = floor((iorbi-1+cellph/2)/cellph) + 1;
binz1 = floor((go+(or*pi/nthet)/2)/(or*pi/nthet)) + 1;
binx2 = binx1 + 1;
biny2 = biny1 + 1;
binz2 = binz1 + 1;
这几句，就是用来计算八个统计区间中心点的坐标的。
在计算前面所讲的统计区间的中心坐标，分配权值之前，我为了处理边缘时程序简洁点，就给那个2*2*9的立体直方图外边又包了一层，形成了一个4*4*11的三维直方图(示意图如下)，原来的2*2*9直方图就是被包在中间的部分。这样，在原来直方图里坐标为(binx=1,biny=2,binz=9)的bin，在新的直方图里坐标为(binx=2,biny=3,binz=10)。
对上面的4*4*11的直方图来个与xoy平面平行的剖面图：
粗实线框就是原三维直方图的剖面，也就是一个block，对于像落在粗实线框与粗虚线框之间的点，其最近邻区间是不够8个的，我为了写程序时省点脑力。。。，就用外扩了的这一圈bin，这样落在粗实线框与粗虚线框之间的统计点有了8个区间，用matlab编程时，那个四层for循环中的部分就只用把那八个累积公式写上，也不用判断是不是在落在像上面粗实线框与粗虚线框之间的那种区域。在程序中2*2*9的直方图为hist3d，4*4*11的直方图为hist3dbig.当在这个hist3dbig中计算都结束后，我把外层这一圈剥去，就是hist3d了。
有了这些准备，我就可以计算出当前像素点的梯度方向幅值应该往hist3dbig中的哪八个bin块累积了。binx1，biny1，binz1 在这里就是那个八个bin块之中离 当前要统计的像素点在直方图中对应的位置 最接近的bin块的下标。binx2，biny2，binz2对应就是最远的bin块的下标了。x1,y1,z1就是bin块(binx1，biny1，binz1)中心点对应的实际像素所在的位置（x1,y1）与梯度方向的角度(z1). 我仍然以原block(即没扩前的block)左上角处作为x1，y1的原点，因为matlab以1作为图像像素索引的开始，我把原点就认为是(1,1)，那(1,1)左边外扩出来的部分，就给以0，-1，-2，-3…这样的坐标,向上也类似，如下图所示，(1,1)位置为红点所示，蓝点处坐标就是(-3,1).
扩展出来的绿块的下标是(binx=1,biny=1,binz1=1),由于像素坐标在红点处为(1,1)，而黄块才是block的第一个cell，对应bin块的下标(2,2).因为下标设计的原因，我在求x1,y1,z1时减了1.5而非0.5.
x1 = (binx1-1.5)*cellpw + 0.5;
y1 = (biny1-1.5)*cellph + 0.5;
z1 = (binz1-1.5)*(or*pi/nthet);
上面的式子中x1，y1还加了0.5，因为像素坐标是离散的，而第一个坐标总是从1开始，这样对如图中第一个cell的中心（黑点）处应该是4.5. z1没加0.5,是因为角度值是从0开始的，并且是连续的。
在signed（即梯度方向从0度到360度）情况下，因为实际上角度的投票区间是首尾相接环形的，若统计间隔是40度，那么0-40度和320-360度就是相邻区间，那么在4*4*11的直方图中，投给binz==11区间(相当于360-380度)的值应该返给binz==2(0-40度)，投给binz==1区间的值应该返给binz==10区间，如4*4*11直方图中所示，对应在程序中就是
if or == 2
hist3dbig(:,:,2) = hist3dbig(:,:,2) + hist3dbig(:,:,nthet+2);
hist3dbig(:,:,(nthet+1)) = hist3dbig(:,:,(nthet+1)) + hist3dbig(:,:,1);
end