HOG算法中用到三线性插值（Trilinear Interpolation）的作用及其用法（包括HOG的matlab源代码分析）

看完第二节，区看后面的代码

1. 写在前面的屁话。

最近的科研任务需要对HOG算法进行魔改，很自然的就需要来看一看算法，看一看代码了。
一开始网上down了一些代码，发现效果很差。分析了一下原因，发现他们都是根据作者论文直接复现的。然后没办法，把matlab中的代码调出来看了一下，果然，事情没有那么简单。

里面分别用了高斯滤波，三线性插值这两个骚操作，简直666.

讲道理，代码我看了两天，真的是没看懂，这也是因为我没有看原始论文的原因，为了图快，也是看了网上的博客，大家都是讲了大概，文章中的一些细节都没讲明白，其中最主要的一点就是，为什么要用三线性插值，怎么用三线性插值。
下面我就把这个做笔记写下来吧。

2. 三线性插值

-----通常将某个变量范围固定划分为几个区域进行某种计算时，由于边界变量与相邻区域也有相关性，如果只对当前区域进行计算而完全忽略与相邻区域的关系，就会产生区域混叠效应。这种混叠效应会在特征向量中产生突变。在这种情况下就需要采用插值算法对计算结果进行修正。在HOG特征提取方法中，位于不同cell交界处的像素如果只对所在的cell进行投影同样会对其他区域产生混叠效应，此时需要采用三维线性插值的方法对梯度向量直方图进行修正。
------插值运算在图像放大中应用较多， 这里借用插值运算的思想对累积直方图进行修正， 即将某个像素点的梯度幅值以不同的权重累加到相应的bin上。确切地说， 这是一种权值分配， 但多数文献和博客中仍然延续 “插值” 这一说法。由于提取HOG特征时需要在两个位置坐标(x, y)和一个方向坐标( θ )共三个维度上进行插值运算， 因此称为三线性插值。
------首先以二维情况为例说明线性插值的思想。在两个维度上进行线性插值称为“双线性插值”。双线性插值在图像放大或图像旋转中应用广泛，它是利用周围4个邻点的灰度值在两个方向上作线性内插以得到待采样点的灰度值，即根据待采样点与相邻点的距离确定相应的权值计算出待采样点的灰度值。如下图所示，点P为待采样点，Q11、Q12、Q21、Q22为点P的四个相邻点，用线性插值法对P点进行插值运算的数学表达式为：

----------------------------------------------------二维线性插值原理图

2.1 在梯度方向上进行线性插值

根据上述插值思想，可在各像素的梯度方向上进行加权运算，如下图所示，将区间[0°,180°]以20°为一个区间划分，每个小区间以中心角度作为直方图的中心数值，假设要对梯度方向为15°的像素点进行处理，显然15与以10和30为中心的直方图最近，应该将该点的梯度幅值加权累加到这两个直方图上，权重分别为(30-15)/20=0.75和(15-10)/20=0.25。

---------------------------------------对待处理像素点的梯度方向进行线性加权示意图
针对这个图，我简单解释一下。在HOG算法中，根据图像梯度方向范围，将其划分为无符号型（0，180°）或者有符号行（-180，180）。我是使用（0，180°）。这样的话，因为默认将角度分为9个方向，因此就事9个bin，每个bin带宽20°，即（0，20），（20，40），…,（160，180） ，因此，在bin中，梯度的每个方向的中心分别是：10，30，50，70，…,170.。

2.2 在像素位置上进行插值

同理运用线性插值方法在各个像素的位置上进行加权运算， 如下图所示， 左图中的方框处为待处理像素点， 它位于block中的C0单元中， 利用该点与四个cell中的中心像素点 （图中4个圆点） 的距离计算权值， 将待处理像素点的梯度幅值分别加权累加到C0、C1、 C2、 C3中相应的直方图上。

----------------------对待处理像素点的位置进行线性加权
这白牛举例分析一下：
------设图中所示像素点的梯度方向是85度，梯度幅值是100，该象素点距离cell中心的左、右、上、下的距离分别为2、6、2、6。首先考虑梯度方向上的插值，若每20度为一个区间，85介于70和90之间，到第四个和第五个角度区间中心的距离分别为15和5，因此若投票值为v，则投票到第四个区间的值是$(5/20)\ast v=0.25v$，投票到第五个区间的值是$(1-1/4)\ast v=0.75v$。接下来考虑在x和y方向上的插值，根据象素点距离各个cell中心的距离，可知在x方向上的权重分配系数为6/8、2/8，在y方向上的权重分配系数也为6/8、2/8。所以梯度幅值分配到第一个cell的值为$100\ast 6/8\ast 6/8=56.25$，分配到第二个cell的值为$100\ast 2/8\ast 6/8=18.75$，分配到第三个cell的值为$100\ast 6/8\ast 2/8=18.75$，分配到第四个cell的值为$100\ast 2/8\ast 2/8=6.25$。
------最后，根据梯度方向上的投票权重，可知：第一个cell的直方图第五个区间得到的投票值为$56.25\ast 0.25=14.0625$，第一个cell的直方图第四个区间得到的投票值 $56.26\ast 0.75=42.1875$；第二个cell的直方图第三个区间得到的投票值为$18.75\ast 0.25=4.6875$，第二个cell的直方图第四个区间得到的投票值为$18.75\ast 0.75=14.0625$；以此类推，可求出第三个和第四个cell的直方图特征。
------上面讲的只是计算block中的一个像素在4个cell的bin中的值，基本都是距离越远，权重越小。

2.3 在像素位置和像素梯度方向上进行三线性插值

综合考虑，在两个位置坐标(x,y)和一个方向坐标(θ)上进行三线性插值，关键要解决的问题是应该在哪些bin上进行加权累加，累加时权值又是多少。将一个像素点处的梯度幅值加权分配到4个cell中与该点梯度方向相邻的2个bin上（这边还需要对梯度幅进行高斯滤波）。按照下面公式修正直方图向量，其中x、y轴表征像素点的空间位置，z轴表征该点的梯度方向（即θ）。对于待处理像素点(x,y)，设其梯度幅值为ω，梯度方向为z，z1和z2分别是与之最邻近的两个bin的中点坐标。梯度直方图h沿x、y、z三个维度的直方图带宽分别为b=[bx,by,bz]，bx=by=8，bz=180°/9。如图6所示为三线性插值计算梯度方向直方图向量的示意图，左图中的方框处为待处理像素点，计算block的每个cell中与该点梯度方向相邻的2个bin，共计8个直方图柱上的权值，将该点的梯度幅值进行加权累加，即形成block中的梯度方向直方图[5]。

由于线性插值法考虑了待采样点周围直接邻点对待采样点的影响， 因此能够克服区域混叠的问题。

--------------------------------------三线性插值计算梯度方向直方图向量

3. 结合代码来分析

为什么要写代码呢，这是因为我们干看上面那8行代码会很懵逼，算法的魅力就在于能推出来公式，还能用代码把公式给实现了。

3.1 关于编程的一些心得

------大一的时候接触了PLC编程，一开始是三菱的，后来是西门子的。讲真，那是土逼的我第一次被自动化的东西震撼到了。没错，就是震撼。
------大二的时候，接触了单片机，这次真的被雷晕了，卧槽，我从，这个更牛逼，比PLC还牛逼，不行，我必须要学会他，然后自己区控制我想控制的东西。于是，从大二开始我便开始学习写代码，那时候学的还是C。当宿舍里的人都在用电脑玩LOL和QQ飞车时，我已经会用单片机驱动四轮小车了。
大四的时候开始接触C++这种让人秃头的东西，然后用这个玩意儿自己搞了一个无人机，四螺旋桨的那种，板子飞控是别人买的，机身时自己焊接的，代码也是自己调试的。（后悔那个时候没去大疆了。）
------研一的时候，2015接触人工智能，当时我跟我硕导都是二把刀子，于是我俩在我硕导的办公室并排坐着，研究了一个月的各种深度学习和经典机器学习的代码，。然后我们总结出来一个结论：如果想看懂一个代码，人后再自己写代码，那么首先样看懂数据的走向。也是就是数据在这个网络模型中时怎么变换的。因此，后来我们实验室后来的师弟师妹们进门的第一步就是学会读取和显示各种类型的数据，然后看懂一个算法的代码，说清楚他们为什么要这样处理数据。

3.2 HOG的主要函数

再该部分，我们主要看的时如下几个函数，因为这涉及HOG处理数据的过程。

function x = normalizeL2Hys(x)
%这个函数是正则化函数，负责对算法中每个像素的梯度幅值进行正则化，这没什么讲的
classToUse = class(x);
x = x./(norm(x,2) + eps(classToUse)); % L2 norm
x(x > 0.2) = 0.2;                     % Clip to 0.2
x = x./(norm(x,2) + eps(classToUse)); % repeat L2 norm
end

function [x1, b1] = computeLowerHistBin(x, binWidth)
%上面我们提到，再构造HOG算法的时候，为了防止窗口滑动的时候出现混叠现象，
%需要在像素位置（x和y轴方向进行插值），同时还要再梯度角度上进行插值，因此需要用到三线性插值。
%这个函数有两个作用：
%1. 对于像素梯度方向的插值，首先我们需要确定其binwidth=20，
%这个函数可以求出来这个方向角度属于哪个bin，以及这个bin的中心角度，通过b1判断这个角度在bin中的位置。
%通过x1判断这个bin的中心角度。
%举个例子，17度，这个梯度方向对应的bin就是1，这个bin的中心角度是15.大家可以带进去看看。
%2. 对于像素位置的插值，首先我们需要确定其binwidth=8，因为我们处理数据是一个一个block来处理的，
%每个block中有4个cells，设置binwidth=8，我们就可以轻松的通过b1判断每个像素属于哪个cell了。
width    = single(binWidth);%转换成单精度
invWidth = 1./width;%求倒数
bin      = floor(x.*invWidth - 0.5);  %floor函数朝负无穷大方向取整，例如-1.2，取值-2。
% min(min(bin))
% max(max(bin))
% Bin 的中心 x1
x1 = width * (bin + 0.5);
%x1(x1==-10) = 10;%weilebaozheng 
% 加2以获得基于1的索引
b1 = int32(bin + 2);
% min(min(b1))
% max(max(b1))
end

function [gMag, gDir] = hogGradient(img,roi)
%这个函数是用于求图像梯度幅值和梯度方向的函数
if nargin == 1
    roi = [];    
    imsize = size(img);
else
    imsize = roi(3:4);
end

img = single(img);

if ndims(img)==3%判断是彩色图像还是灰度图像。三维表示彩色图像
    rgbMag = zeros([imsize(1:2) 3], 'like', img);
    rgbDir = zeros([imsize(1:2) 3], 'like', img);
    
    for i = 1:3
        [rgbMag(:,:,i), rgbDir(:,:,i)] = computeGradient(img(:,:,i),roi);
    end
    
    % find max color gradient for each pixel
    [gMag, maxChannelIdx] = max(rgbMag,[],3);
    
    % extract gradient directions from locations with maximum magnitude
    sz = size(rgbMag);
    [rIdx, cIdx] = ndgrid(1:sz(1), 1:sz(2));
    ind  = sub2ind(sz, rIdx(:), cIdx(:), maxChannelIdx(:));
    gDir = reshape(rgbDir(ind), sz(1:2));
else%如果是灰度图像则是直接调用这个函数，下面我来介绍这个函数
    [gMag,gDir] = computeGradient(img,roi);
end

end

function [gMag,gDir] = computeGradient(img,roi)

if isempty(roi)%首先判断是不是要求感兴趣区域的图像像素，如果没有输入图像区域，则表示跳过这步。
%这些代码我就不讲了，关于求梯度很简单的。
    gx = zeros(size(img), 'like', img);
    gy = zeros(size(img), 'like', img);
    
    gx(:,2:end-1) = conv2(img, [1 0 -1], 'valid');
    gy(2:end-1,:) = conv2(img, [1;0;-1], 'valid');
    
    % forward difference on borders
    gx(:,1)   = img(:,2)   - img(:,1);
    gx(:,end) = img(:,end) - img(:,end-1);
    
    gy(1,:)   = img(2,:)   - img(1,:);
    gy(end,:) = img(end,:) - img(end-1,:);
else%如果是求整个图像的所有像素的梯度，则直接调用这个函数，下面我们再看看这个函数
    [gx, gy] = computeGradientROI(img, roi);
end

% return magnitude and direction
gMag = hypot(gx,gy);
gDir = atan2d(-gy,gx);
end

function [gx, gy] = computeGradientROI(img, roi)
img    = single(img);
imsize = size(img);

% roi is [r c height width]
rIdx = roi(1):roi(1)+roi(3)-1;
cIdx = roi(2):roi(2)+roi(4)-1;

imgX = coder.nullcopy(zeros([roi(3)   roi(4)+2], 'like', img)); %#ok<NASGU>
imgY = coder.nullcopy(zeros([roi(3)+2 roi(4)  ], 'like', img)); %#ok<NASGU>

% replicate border pixels if ROI is on the image border. 
if rIdx(1) == 1 || cIdx(1)==1  || rIdx(end) == imsize(1) ...
        || cIdx(end) == imsize(2)
    
    if rIdx(1) == 1
        padTop = img(rIdx(1), cIdx);
    else
        padTop = img(rIdx(1)-1, cIdx);
    end
    
    if rIdx(end) == imsize(1)
        padBottom = img(rIdx(end), cIdx);
    else
        padBottom = img(rIdx(end)+1, cIdx);
    end
    
    if cIdx(1) == 1
        padLeft = img(rIdx, cIdx(1));
    else
        padLeft = img(rIdx, cIdx(1)-1);
    end
    
    if cIdx(end) == imsize(2)
        padRight = img(rIdx, cIdx(end));
    else
        padRight = img(rIdx, cIdx(end)+1);
    end
    
    imgX = [padLeft img(rIdx,cIdx) padRight];
    imgY = [padTop; img(rIdx,cIdx);padBottom];
else  
    imgX = img(rIdx,[cIdx(1)-1 cIdx cIdx(end)+1]);
    imgY = img([rIdx(1)-1 rIdx rIdx(end)+1],cIdx);
end

gx = conv2(imgX, [1 0 -1], 'valid');
gy = conv2(imgY, [1;0;-1], 'valid');
end

好了，下面开始这长的主要环节了，即分析算法是如何提取HOG特征，以及数据再HOG算法中的变换方式。

闲话少数，直接上代码：

function [cellhog,hog] = extractHOG(gMag, gDir, gaussianWeights, weights, params)
    %----------------------------------------------------------------------
    %原始函数是没有cellhog这个变量的，因为需要用到这个变量，因此把他提取出来了。
    %原始的算法只输出hog这一个变量，并且把他reshape成了一列向量。
    %这里输入的gMag和gDir都是与图像img相同尺寸的大小，
    %gMag是像素的梯度幅值，gDir是像素的梯度方向。
    %里面每个数据分别是该图像每个像素的梯度幅值和梯度方向。因此是个矩阵
    %
    %----------------------------------------------------------------------
if isempty(coder.target)         
    hog = visionExtractHOGFeatures(gMag, gDir, gaussianWeights, params, weights);                  
else
    featureClass = 'single';%定义特征类型为单精度
    
    if params.UseSignedOrientation%默认为0，因为我们习惯了使用180度
        % 定义gDir的范围在 [0 360]之间
        histRange = single(360);%单精度型变量
    else
        % 把方向转换成无符号型，范围在 [0 180]之间
        histRange = single(180);%单精度型变量，这个histRange是像素梯度方向的的范围
    end
    
    %如果 gDir 范围是 [-180 180], 转换到  [0 180] 或者 [0 360]之间
    negDir = gDir < 0;%真为1，假为0
    gDir(negDir) = histRange + gDir(negDir);
    
    % 为所有的cells定位方向bin
    binWidth = histRange/cast(params.NumBins, featureClass);
    %例如，计算在180度的情况下，有9个bin，即NumBins=9，则每个bin的宽度是20度。
    %这里的NumBins可以根据需要自行设定，但是HOG的原文说这9好使那就是9吧。
    [x1, b1] = computeLowerHistBin(gDir, binWidth);
    %利用该公式计算出当前像素的梯度方向在bin中哪个位置。x1表示当前梯度方向角度整数角度中间值，
    %比如第一个bin是0-20度，那么这个值就是10，b1表示该角度在bin中的第几格。
    wDir = 1 - (gDir - x1)./binWidth;
    %方向权重，与中间值偏离越大，权重越小。这是很自然的。
    %不理解的话可以自己带入参数试试。这是一个与图像大小相同的矩阵。binWidth=20
    %这边其实就是对像素的梯度方向进行线性插值********记住哦，这边已经开始处理数据了。
    
    blockSizeInPixels = params.CellSize.*params.BlockSize;
    %（8，8）*（2，2），表示block包含像素的数量（16，16），
    blockStepInPixels = params.CellSize.*(params.BlockSize - params.BlockOverlap);
    %表示block每次移动几个像素。
    
    r = 1:blockSizeInPixels(1);%行1，2，3，4，....16
    c = 1:blockSizeInPixels(2);%列1，2，3，4，....16
    
    nCells  = params.BlockSize;%每个block的尺寸，用于计算每个block中cell的个数是 2*2
    nBlocks = vision.internal.hog.getNumBlocksPerWindow(params);%每个窗口中block的数量
    numCellsPerBlock = nCells(1)*nCells(2);%每个block中cell的个数,4
    
    %下面就是重点了，首先，作者定义了一个变量hog，用于存贮最终的hog特征，让我们来分析一下他的思想
    hog = coder.nullcopy(zeros([params.NumBins*numCellsPerBlock, nBlocks],featureClass)); 
    %（bin的数量乘以每个block中cell的数量，nBlocks是窗口中block的数量）。
    %hog的存储思想是：每列存储一个block中4个cell的bin，每个bin占9个单位，每个cell一个bin。
    %（9*4,n），n是每个窗口中block的个数。每个cell中有9个bin，每个block中有4个cell，所以每个block中有36个bin空间，分别是4个cell的bin
    %nBlocks = 15*7
    %定义的这个hog，其思想就是，每一列就是一个block所包含的bin的数量，1到9是cell（1，1）的bin，
    %10到18是cell(1,2)的bin，19到27是cell(2,1)的bin，28到36是cell(2,2)的bin。
    % 扫描窗口内的所有block，起到遍历图像的目的
    %先来说一下遍历图像的思想：外层循环控制行，内层循环控制列。
    %由此可以看出，算法中的block是先横向滑动，然后再换到第二行再横向滑动，从而遍历整张图像。
    for j = 1:nBlocks(2);       
        for i = 1:nBlocks(1)            
            %下面是对方向权重进行线性插值加权
            wz1 = wDir(r,c); %wDir是对梯度方向进行线性加权后的梯度方向值。
            %算法通过对r和c进行递加运算获取每个block对应的图像数据。
            %wz1 这是一个与block大小相同的矩阵
            w = trilinearWeights(wz1, weights); %以这边需要去看下面的两个函数。这一步就是实现了对三个变量的三线性加权                       
            % 对幅值进行高斯滤波加权
            m = gMag(r,c) .* gaussianWeights; %对像素幅值进行高斯平滑处理。这是一个与block大小相同的矩阵           
            % interpolate magnitudes for binning 对梯度幅值进行插值并存到bin中
            
            %下面的步骤完成了将梯度幅值和梯度方向进行加权融合的任务，从而将会两个数据整合成一个数据，
            mx1y1z1 = m .* w.x1_y1_z1;
            mx1y1z2 = m .* w.x1_y1_z2;
            mx1y2z1 = m .* w.x1_y2_z1;
            mx1y2z2 = m .* w.x1_y2_z2;
            mx2y1z1 = m .* w.x2_y1_z1;
            mx2y1z2 = m .* w.x2_y1_z2;
            mx2y2z1 = m .* w.x2_y2_z1;
            mx2y2z2 = m .* w.x2_y2_z2;            
            %--------------------------------------------------------------
            %下面开始处理该存储在哪个bin里面
            %--------------------------------------------------------------
            orientationBins = b1(r,c);%获取该block区域中像素的梯度方向再对应bin中的位置。         
            % 将块直方图进行初始化为零。
            %注意了，这又是一个关键变量h，h适用于存储每个block中特想的hog特征值的。
            %也就是每循环一次，h都会把它里面的值存储到hog中，到循环最后，再把hog  reshape成一个一维向量
            h = zeros(params.NumBins+2, nCells(1)+2, nCells(2)+2, featureClass);%单精度，h用于存储每个block中的bin。
            %再默认的情况下：
            %h = zeros(11, 4, 4, featureClass);
            %--------------------------------------------------------------
            %这边需要着重搞懂直方图的存储策略
            %h被定义成一个由4个11*4的三维矩阵，数据类型为单精度类型
            %--------------------------------------------------------------            
            % 将内插量值放到到块直方图中,横着计算block中的每个像素。将每个像素的梯度方向加权后的结果放进bin中。
            for x = 1:blockSizeInPixels(2)%1：16,用列数，表示横着走，
                cx = weights.cellX(x);%判断这个像素属于哪个cell
                for y = 1:blockSizeInPixels(1)%1：16
                    z  = orientationBins(y,x);
                    %这个像素的梯度方向属于哪个bin，这里判断出来的角度事实际属于的bin，而z+1是比这个z大且相邻的另一个bin。
                    %因为算法需要在两个相邻bin内进行加权，根据距离进行加权。
                    cy = weights.cellY(y);%1，2，3
                    %针对这边的解释，我们可以参考上面第二节中的那一大坨公示了，是不是一样的。
                    %单纯这样看，好多人可能还不是很明白，为什么会这样写。我们给h里面的每个变量赋值，大家就能看明白了。
                    %在这里，
                    h(z,   cy,   cx  ) = h(z,   cy,   cx  ) + mx1y1z1(y,x);%这边其实是一个实数，不是三维矩阵
                    h(z+1, cy,   cx  ) = h(z+1, cy,   cx  ) + mx1y1z2(y,x);
                    h(z,   cy+1, cx  ) = h(z,   cy+1, cx  ) + mx1y2z1(y,x);
                    h(z+1, cy+1, cx  ) = h(z+1, cy+1, cx  ) + mx1y2z2(y,x);
                    h(z,   cy,   cx+1) = h(z,   cy,   cx+1) + mx2y1z1(y,x);
                    h(z+1, cy,   cx+1) = h(z+1, cy,   cx+1) + mx2y1z2(y,x);
                    h(z,   cy+1, cx+1) = h(z,   cy+1, cx+1) + mx2y2z1(y,x);
                    h(z+1, cy+1, cx+1) = h(z+1, cy+1, cx+1) + mx2y2z2(y,x);
                end
            end
            
            % h装填，这里这样操作是因为上面创建h时在x和y上+2，之所以加2，是因为我们在计算像素属于哪个cell时，出现1，2，3，然而3不是我们想要的。
            h(2,:,:)     = h(2,:,:)     + h(end,:,:);
            h(end-1,:,:) = h(end-1,:,:) + h(1,:,:);
            
            % 仅保留块直方图的有效部分
            %这里解释一下，为什么其他部分是无效的直方图。首先，
            h = h(2:end-1,2:end-1,2:end-1);
            
            %正则化并将块添加到特征向量         
            hog(:,i,j) = normalizeL2Hys(h(:));%从代码上观察，每个h表示一个block中的所有bin，在hog中，每个j表示一页，每一页的矩阵便是一个block的所有bin值，       
            r = r + blockStepInPixels(1);
        end
        r = 1:blockSizeInPixels(1);
        c = c + blockStepInPixels(2);
    end
    cellhog = hog;
    hog = reshape(hog, 1, []);%把hog的特征reshape一个一行N列的向量，用于到SVM等算法中训练
end

4. 参考文献

特此感谢
黄冬丽, 戴健文, 冯超, et al. HOG特征提取中的三线性插值算法[J]. 电脑知识与技术, 2012(31):7548-7551.