HOG特征matlab代码实现

%Matlab版HOG代码
function F = hogcalculator（img cellpw cellph nblockw nblockh...
    nthet overlapissigned normmethod）;

% HOG特征由Dalal在2005 cvpr 的一篇论文中提出

% NORMMETHOD：重叠块中的特征标准化函数的方法
%       e为一个设定的很小的数使分母不为0
%       v为标准化前的特征向量
%       ‘none‘ which means non-normalization;
%       ‘l1‘ which means L1-norm normalization; V=V/（V+e）
%       ‘l2‘ which means L2-norm normalization; V=V/根号（V平方+e平方）
%       ‘l1sqrt‘V=根号（V/（V+e））
%       ‘l2hys‘l2的省略形式。将V最大值限制为0.2

if nargin < 2
    % 在DALAL论文中指出的在rows:128*columns:64情况下的最佳值，设定为DEFAULT
    cellpw = 8;
    cellph = 8;
    nblockw = 2;
    nblockh = 2;
    nthet = 9;
    overlap = 0.5;
    issigned = ‘unsigned‘;
    normmethod = ‘l2hys‘;
else
    if nargin < 9
        error（‘输入参数不足.‘）;
    end
end

[M N K] = size（img）;  %M为行数，N为列数，K为维数
if mod（Mcellph*nblockh） ~= 0   %行数必须为块的高度的整数倍
    error（‘图片行数必须为块的高度的整数倍.‘）;
end
if mod（Ncellpw*nblockw） ~= 0   %列数必须为块的宽度的整数倍
    error（‘图片列数必须为块的宽度的整数倍.‘）;
end                             
if mod（（1-overlap）*cellpw*nblockw cellpw） ~= 0 ||...  %要使滑步后左边是整数
        mod（（1-overlap）*cellph*nblockh cellph） ~= 0
    error（‘滑步的像素个数必须为细胞单元尺寸的整数倍‘）;
end

%设置高斯空间权值窗口的方差
delta = cellpw*nblockw * 0.5;


%计算梯度矩阵  梯度的计算【-1，0，1】效果是很好的，而3*3的sobel算子或者2*2的对角矩阵反而会系统的降低效果
hx = [-101];
hy = -hx‘;   %转置
gradscalx = imfilter（double（img）hx）;  %imfilter是滤波器，hx表示滤波掩膜
gradscaly = imfilter（double（img）hy）;

if K > 1
    gradscalx = max（max（gradscalx（::1）gradscalx（::2）） gradscalx（::3））;  %取RGB中最大值
    gradscaly = max（max（gradscaly（::1）gradscaly（::2）） gradscaly（::3））;
end
gradscal = sqrt（double（gradscalx.*gradscalx + gradscaly.*gradscaly））;  %梯度矩阵 gradscal

% 计算梯度方向矩阵
gradscalxplus = gradscalx+ones（size（gradscalx））*0.0001;  %防止为0，所以gradscalx加了0.0001
gradorient = zeros（MN）;                                 %初始化梯度方向矩阵
% unsigned situation: orientation region is 0 to pi.
if strcmp（issigned ‘unsigned‘） == 1                     %无向的情况
    gradorient =...
        atan（gradscaly./gradscalxplus） + pi/2;           %加pi/2因为atan的区间取值从-pi/2开始
    or = 1;
else
    % signed situation: orientation region is 0 to 2*pi. %有向的情况
    if strcmp（issigned ‘signed‘） == 1
        idx = find（gradscalx >= 0 & gradscaly >= 0）;
        gradorient（idx） = atan（gradscaly（idx）./gradscalxplus（idx））;
        idx = find（gradscalx < 0）;
        gradorient（idx） = atan（gradscaly（idx）./gradscalxplus（idx）） + pi;
        idx = find（gradscalx >= 0 & gradscaly < 0）;
        gradorient（idx） = atan（gradscaly（idx）./gradscalxplus（idx）） + 2*pi;
        or = 2;
    else
     %  error（‘Incorrect ISSIGNED parameter.‘）;
        error（‘参数ISSIGNED输入有误‘）;
    end
end

% 计算块的滑步
xbstride = cellpw*nblockw*（1-overlap）;   %x方向的滑步
ybstride = cellph*nblockh*（1-overlap）;
xbstridend = N - cellpw*nblockw + 1;     %x方向块左角能达到的最大值
ybstridend = M - cellph*nblockh + 1;

% 块总数=ntotalbh*ntotalbw
ntotalbh = （（M-cellph*nblockh）/ybstride）+1; %除了第一个，后面每