Matlab 提取Hog特征方法详细解析

     前段时间，写了一篇《matlab Hog特征提取浅析》，当时只是完成了对Hog特征提取方法的简单原理分析。最近，我深入地研究了Hog特征的提取方法，并运用积分直方图对其加速。为了便于以后运用，特将这一过程中的构思，方法，问题一一记录下来。

   相关资料：

   1.http://blog.csdn.net/zouxy09/article/details/7929348/
   2.http://blog.csdn.net/xiaowei_cqu/article/details/17928733
   3.http://blog.csdn.net/WLMLBY/article/details/53446943

   提取Hog特征的方法是先把窗口图像分割成cell，再把cell分割成block。然后提取block中各像素点的方向梯度合成直方图，最后把这些block的直方图组合起来就得到了图像的Hog特征。
   下面具体说一下我的编程思路。

   主程序(Main.m)：

close all;clear all;clc;filePath = 'E:\HOG\HogSamples\';fileExt = '*.bmp';files = dir(fullfile(filePath,fileExt));len = size(files,1);% Generate a test Image.% testImg = reshape((1:16:256),4,4);% hogFeature = ExtrHogFeature(testImg,3,3,1,180,9);for i = 1:len    tic;    disp(['----------',num2str(i),'----------'])    srcImg = imread(strcat(filePath,files(i,1).name));    hogFeature = ExtrHogFeature(srcImg,5,4,2,360,15);    toc;end

   在主程序中，我遍历了指定文件夹下(“E:\HOG\HogSamples\”)的指定后缀名（*.bmp）的文件。并将遍历得到的图像作为窗口图像传递给了ExtrHogFeature函数。ExtrHogFeature函数就是我编写的，用来提取Hog特征的函数。    ExtrHogFeature函数：

function feature = ExtrHogFeature(imgFile,cellStep,blockStep,overlap,angle,binNum)if nargin < 2    cellStep = 3;    blockStep = 4;    overlap = 1;    angle = 180;    binNum = 9;elseif nargin < 6    error('Input parameters are not enough.');elseif overlap < 0 || overlap >= blockStep    error('The input parameter "overlap" is wrong!');elseif angle ~= 180 && angle ~= 360    error('The input parameter "angle" is only 180 or 360!');end[h,w,d] = size(imgFile);if(d ~= 1)    grayImg = rgb2gray(imgFile);else    grayImg = imgFile;end% Correct the width and height of imgFile.if w < cellStep    cellNum_W = 1;else    cellNum_W = ceil(w / cellStep);    if cellNum_W < blockStep        blockNum_W = 1;    else        blockNum_W = floor((cellNum_W - blockStep) / (blockStep - overlap)) + 1;    endendcorrW = cellStep * ((blockStep - overlap) * (blockNum_W - 1) + blockStep);if h < cellStep    cellNum_H = 1;else    cellNum_H = ceil(h / cellStep);    if cellNum_H < blockStep        blockNum_H = 1;    else        blockNum_H = floor((cellNum_H - blockStep) / (blockStep - overlap)) + 1;    endendcorrH = cellStep * ((blockStep - overlap) * (blockNum_H - 1) + blockStep);% Bilinear interpolationcorrImg = BiliInter(grayImg,corrH,corrW);% The following contents need to be modified.corrImg = double(corrImg);flipImg = FlipImg(corrImg,1,1);Gx = zeros(corrH,corrW);Gy = zeros(corrH,corrW);range = angle / binNum;for i = 1:corrH    for j = 1:corrW        Gx(i,j) = flipImg(i,j + 1 + 1) - flipImg(i,j);        Gy(i,j) = flipImg(i + 1 + 1,j) - flipImg(i,j);        Gxy(i,j) = sqrt((Gx(i,j))^2 + (Gy(i,j))^2);    endendif angle == 180    Axy = (atan(Gy ./ Gx) / pi) * 180;    Axy(find(Axy < 0)) = Axy(find(Axy < 0)) + 180;elseif angle == 360    Axy = (atan2(Gx,Gy) / pi) * 180;    Axy(find(Axy < 0)) = Axy(find(Axy < 0)) + 360;endAxy(isnan(Axy)) = 0;Axy = floor((Axy / range) + 1);coluHist = zeros(1,corrW,binNum);inteHist = zeros(corrH,corrW,binNum);for i = 1:corrW    coluHist(1,i,Axy(1,i)) = Gxy(1,i);    inteHist(1,i,Axy(1,i)) = Gxy(1,i);    if i > 1        inteHist(1,i,:) = inteHist(1,i,:) + inteHist(1,i - 1,:);    endendfor i = 2:corrH    for j = 1:corrWcoluHist(1,j,Axy(i,j)) = coluHist(1,j,Axy(i,j)) + Gxy(i,j);        if j == 1                      inteHist(i,j,:) = coluHist(1,j,:);        else                      inteHist(i,j,:) = inteHist(i,j - 1,:) + coluHist(1,j,:);        end    endendtmpFeature = zeros(1,binNum);for i = 1:blockNum_H    for j = 1:blockNum_W        i1 = (i * (blockStep - overlap) + overlap - blockStep + 1) * cellStep;        j1 = (j * (blockStep - overlap) + overlap - blockStep + 1) * cellStep;        i2 = (i * (blockStep - overlap) + overlap - blockStep + 1) * cellStep;        j2 = (j * (blockStep - overlap) + overlap) * cellStep;        i3 = (i * (blockStep - overlap) + overlap) * cellStep;        j3 = (j * (blockStep - overlap) + overlap - blockStep + 1) * cellStep;        i4 = (i * (blockStep - overlap) + overlap) * cellStep;        j4 = (j * (blockStep - overlap) + overlap) * cellStep;        if i == 1 && j == 1            tmpFeature(:) = inteHist(i4,j4,:);            feature = tmpFeature;        else            tmpFeature(:) = inteHist(i4,j4,:) - inteHist(i2,j2,:) - inteHist(i3,j3,:) + inteHist(i1,j1,:);            feature = [feature,tmpFeature];        end    endendmaxValue = max(feature(:));minValue = min(feature(:));feature = (feature - minValue) / (maxValue - minValue) * 1000;% save ExtrHogFeature_Data

   在ExtrHogFeature函数中：

   1.参数意义：

   (1).imgFile：输入的窗口图像文件。

   (2).cellStep：cell的步长，即一个cell有多少个像素。理论上，cell的宽和高是可以不一样的，但是在实际编写程序时，这么做没有多大的意义，单纯增加代码量。所以，我只设置一个cellStep，默认cell是正方形。同理，下面的blockStep和overlap也是这样。

   (3).blockStep：block的步长，即一个block有多少个cell。注意：不是有多少个像素，是有多少个cell。

   (4).overlap：重叠数目。即后（下）一个block与前（上）一个block有几列(行)cell重叠。

   (5).angle：只能是两个值，180或者360。

   (6).binNum：积分直方图的区域数量。

    2.判断输入参数是否正确。

   (1).判断参数的数量。参数的数量只能是两种，1个或者6个。如果是1个，对除imgFile之外的其他参数进行默认赋值。

   (2).判断overlap是否正确。因为overlap表示两个block之间的重叠数量。所以，overlap不能大于或者等于blockStep。

   (3).判断angle是否正确。angle只能是180或360。

    3.判断图像的维数，并进行灰度化。这一步对应着相关资料1中的“标准化gamma空间和颜色空间”。从数学原理上看，只是对图像矩阵进行了开方。为了简化计算，对窗口图像进行灰度化。

    4.校正窗口图像的宽和高。在提取图像特征时，我一直秉持“不能减少图像的特征”的原则，即不能对图像进行裁剪。因为原始图像的宽和高，不一定能生成整数量的cell，更不一定能生成整数量的block。所以，需要利用双线性插值对原始图像进行校正（即缩放），使能生成整数量的block。

    举个例子：ExtrHogFeature(testImg,3,4,2,180,9),其中testImg是宽和高均为20个像素的测试窗口图像。cellStep=3，cellNum_W（每行cell的个数）=20/3=6.6，取cellNum_W =ceil(20/3)=7。

    设blockNum_W是每一行block的数量，满足下列公式：

    (blockStep - overlap) * (blockNum_W - 1) + blockStep >= cellNum_W   ①

    取满足公式的最小整数。即：

    blockNum_W = floor((cellNum_W - blockStep) / (blockStep - overlap)) + 1   ②

    由①式计算得出，blockNum_W>=2.5，所以（也可由②式得）blockNum_W=3。cellNum_W=（4-2）*（3-1）+4=8。所以，经过校正后的图像corrImg的宽corrW=8*3=24个像素。corrH计算同理。

    根据corrW和corrH，对原始图像进行双线性插值的缩放。双线性插值的函数BiliInter在本篇的末尾。因为不是本篇的主要内容，且易于理解，所以不做额外的讲解。

    5.计算每个像素的方向梯度。水平模板fx=[-1,0,1]，垂直模板fy=fx’。

为了计算corrImg的4条边上的方向梯度，要将4条边翻转。同双线性插值一样，翻转图像的函数FlipImg也在本篇的末尾。

    注意：atan和atan2的区别。atan（y/x）根据正切值y/x求出对应的角度，是2象限的反正切。atan2（y,x）的取值不仅取决于y/x的正切值，还取决于(x,y)落于哪个象限，是4象限的反正切。所以，当angle=180的时候，经过atan计算小于0的值是在-pi/2~0之间，需要+pi，把反正切值校正到pi/2~pi。当angle=360的时候，经过atan2计算小于0的值是在-pi~0之间，需要+2pi，把反正切值校正到pi~2pi。

    Axy中可能会出现0/0=Nan的情况,造成程序中止。所以，要有Axy(isnan(Axy)) = 0。

    6.计算原始窗口图像的积分直方图inteHist。我是根据相关资料2，编写整理的代码。将资料中的2维的积分图计算方法做了适当的修改，应用到3维的积分直方图计算。

    7.利用积分直方图，计算窗口图像的每个block的积分直方图，把这些block的直方图组合起来就得到了图像的Hog特征。并归一化到0~1000的范围内。

   BiliInter函数：

function img = BiliInter(imgFile,dstH,dstW)if dstH <= 0 || dstW <= 0    error('Input parameters are wrong!');end[srcH,srcW] = size(imgFile);img = zeros(dstH,dstW);for i = 0:dstH - 1    for j = 0:dstW - 1        actualW = (j / (dstW - 1)) * (srcW - 1);        actualH = (i / (dstH - 1)) * (srcH - 1);        w_Scale = actualW - floor(actualW);        h_Scale = actualH - floor(actualH);        leftUp = [floor(actualH) + 1,floor(actualW) + 1];        rightUp = [floor(actualH) + 1,ceil(actualW) + 1];        leftLow = [ceil(actualH) + 1,floor(actualW) + 1];        rightLow = [ceil(actualH) + 1,ceil(actualW) + 1];        tmpUp = (imgFile(rightUp(1),rightUp(2)) - imgFile(leftUp(1),leftUp(2))) * w_Scale + imgFile(leftUp(1),leftUp(2));        tmpLow = (imgFile(rightLow(1),rightLow(2)) - imgFile(leftLow(1),leftLow(2))) * w_Scale + imgFile(leftLow(1),leftLow(2));        img(i + 1,j + 1) = (tmpLow - tmpUp) * h_Scale + tmpUp;    endendimg = uint8(img);

   FlipImg函数：

function img = FlipImg(imgFile,flipRowNum,flipColNum)% The flipRowNum is added on imgFile one way.% The flipColNum is similar to flipRowNum.[srcH,srcW] = size(imgFile);if flipRowNum > srcH || flipColNum > srcW    error('Input parameters are too big!');enddstH = srcH + 2 * flipRowNum;dstW = srcW + 2 * flipColNum;tmpImgRow = zeros(dstH,srcW);img = zeros(dstH,dstW);% Fliped with Row.tmpImgRow(1:flipRowNum,:) = imgFile(flipRowNum:-1:1,:);tmpImgRow(flipRowNum + 1:flipRowNum + srcH,:) = imgFile(:,:);tmpImgRow(flipRowNum + srcH + 1:end,:) = imgFile(srcH:-1:srcH - flipRowNum + 1,:);% Fliped with Column.img(:,1:flipColNum) = tmpImgRow(:,flipColNum:-1:1);img(:,flipColNum + 1:flipColNum + srcW) = tmpImgRow(:,:);img(:,flipColNum + srcW + 1:end) = tmpImgRow(:,srcW:-1:srcW - flipColNum + 1);