function F = hogcalculator(img, cellpw, cellph, nblockw, nblockh,...
nthet, overlap, isglobalinterpolate, issigned, normmethod)
% HOGCALCULATOR calculate R-HOG feature vector of an input image using the
% procedure presented in Dalal and Triggs's paper in CVPR 2005.
%
% Author: timeHandle
% Time: March 24, 2010
% May 12，2010 update.
%
% this copy of code is written for my personal interest, which is an
% original and inornate realization of [Dalal CVPR2005]'s algorithm
% without any optimization. I just want to check whether I understand
% the algorithm really or not, and also do some practices for knowing
% matlab programming more well because I could be called as 'novice'.
% OpenCV 2.0 has realized Dalal's HOG algorithm which runs faster
% than mine without any doubt, ╮(╯▽╰)╭ . Ronan pointed a error in
% the code，thanks for his correction. Note that at the end of this
% code, there are some demonstration code，please remove in your work.
%
% F = hogcalculator(img, cellpw, cellph, nblockw, nblockh,
% nthet, overlap, isglobalinterpolate, issigned, normmethod)
%
% IMG:
% IMG is the input image.
%
% CELLPW, CELLPH:
% CELLPW and CELLPH are cell's pixel width and height respectively.
%
% NBLOCKW, NBLCOKH:
% NBLOCKW and NBLCOKH are block size counted by cells number in x and
% y directions respectively.
%
% NTHET, ISSIGNED:
% NTHET is the number of the bins of the histogram of oriented
% gradient. The histogram of oriented gradient ranges from 0 to pi in
% 'unsigned' condition while to 2*pi in 'signed' condition, which can
% be specified through setting the value of the variable ISSIGNED by
% the string 'unsigned' or 'signed'.
%
% OVERLAP:
% OVERLAP is the overlap proportion of two neighboring block.
%
% ISGLOBALINTERPOLATE:
% ISGLOBALINTERPOLATE specifies whether the trilinear interpolation
% is done in a single global 3d histogram of the whole detecting
% window by the string 'globalinterpolate' or in each local 3d
% histogram corresponding to respective blocks by the string
% 'localinterpolate' which is in strict accordance with the procedure
% proposed in Dalal's paper. Interpolating in the whole detecting
% window requires the block's sliding step to be an integral multiple
% of cell's width and height because the histogram is fixing before
% interpolate. In fact here the so called 'global interpolation' is
% a notation given by myself. at first the spatial interpolation is
% done without any relevant to block's slide position, but when I was
% doing calculation while OVERLAP is 0.75, something occurred and
% confused me o__O"… . This let me find that the operation I firstly
% did is different from which mentioned in Dalal's paper. But this
% does not mean it is incorrect ^◎^, so I reserve this. As for name,
% besides 'global interpolate', any others would be all ok, like 'Lady GaGa'
% or what else, :-).
%
% NORMMETHOD：
% NORMMETHOD is the block histogram normalized method which can be
% set as one of the following strings:
% 'none', which means non-normalization;
% 'l1', which means L1-norm normalization;
% 'l2', which means L2-norm normalization;
% 'l1sqrt', which means L1-sqrt-norm normalization;
% 'l2hys', which means L2-hys-norm normalization.
% F：
% F is a row vector storing the final histogram of all of the blocks
% one by one in a top-left to bottom-right image scan manner, the
% cells histogram are stored in the same manner in each block's
% section of F.
%
% Note that CELLPW*NBLOCKW and CELLPH*NBLOCKH should be equal to IMG's
% width and height respectively.
%
% Here is a demonstration, which all of parameters are set as the
% best value mentioned in Dalal's paper when the window detected is 128*64
% size(128 rows, 64 columns):
% F = hogcalculator(window, 8, 8, 2, 2, 9, 0.5,
% 'localinterpolate', 'unsigned', 'l2hys');
% Also the function can be called like:
% F = hogcalculator(window);
% the other parameters are all set by using the above-mentioned "dalal's
% best value" as default.
%
if nargin < 2
% set default parameters value.
cellpw = 8;
cellph = 8; % 原来是8x8!!!!!!!!!!!!!!!!!!!!!!!!!!
nblockw = 2;
nblockh = 2;
nthet = 9;
overlap = 0.5;
isglobalinterpolate = 'localinterpolate';
issigned = 'unsigned';
normmethod = 'l2hys';
else
if nargin < 10
error('Input parameters are not enough.');
end
end
% check parameters's validity.
[M, N, K] = size(img);
if mod(M,cellph*nblockh) ~= 0
error('IMG''s height should be an integral multiple of CELLPH*NBLOCKH.');
end
if mod(N,cellpw*nblockw) ~= 0
error('IMG''s width should be an integral multiple of CELLPW*NBLOCKW.');
end
if mod((1-overlap)*cellpw*nblockw, cellpw) ~= 0 ||...
mod((1-overlap)*cellph*nblockh, cellph) ~= 0
str1 = 'Incorrect OVERLAP or ISGLOBALINTERPOLATE parameter';
str2 = ', slide step should be an intergral multiple of cell size';
error([str1, str2]);
end
% set the standard deviation of gaussian spatial weight window.
delta = cellpw*nblockw * 0.5;
% calculate gradient scale matrix.
hx = [-1,0,1];
hy = -hx';
gradscalx = imfilter(double(img),hx);
gradscaly = imfilter(double(img),hy);
%
if K > 1
maxgrad = sqrt(double(gradscalx.*gradscalx + gradscaly.*gradscaly));
[gradscal, gidx] = max(maxgrad,[],3);
gxtemp = zeros(M,N);
gytemp = gxtemp;
for kn = 1:K
ttempx = gradscalx(:,:,kn);
ttempy = gradscaly(:,:,kn);
tmpindex = find(gidx==kn);
gxtemp(tmpindex) = ttempx(tmpindex);
gytemp(tmpindex) =ttempy(tmpindex);
end
gradscalx = gxtemp;
gradscaly = gytemp;
else
gradscal = sqrt(double(gradscalx.*gradscalx + gradscaly.*gradscaly));
end
% calculate gradient orientation matrix.
% plus small number for avoiding dividing zero.
gradscalxplus = gradscalx+ones(size(gradscalx))*0.0001;
gradorient = zeros(M,N);
% unsigned situation: orientation region is 0 to pi.
if strcmp(issigned, 'unsigned') == 1
gradorient =...
atan(gradscaly./gradscalxplus);
gradorient(gradorient<0) = gradorient(gradorient<0)+pi;
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.');
end
end
% calculate block slide step.
xbstride = cellpw*nblockw*(1-overlap);
ybstride = cellph*nblockh*(1-overlap);
xbstridend = N - cellpw*nblockw + 1;
ybstridend = M - cellph*nblockh + 1;
% calculate the total blocks number in the window detected, which is
% ntotalbh*ntotalbw.
ntotalbh = ((M-cellph*nblockh)/ybstride)+1;
ntotalbw = ((N-cellpw*nblockw)/xbstride)+1;
% generate the matrix hist3dbig for storing the 3-dimensions histogram. the
% matrix covers the whole image in the 'globalinterpolate' condition or
% covers the local block in the 'localinterpolate' condition. The matrix is
% bigger than the area where it covers by adding additional elements
% (corresponding to the cells) to the surround for calculation convenience.
if strcmp(isglobalinter
