matlab intergral,matlab學習：人臉識別之HOG（Histograms of Oriented Gradients）
HOG descriptors 是應用在計算機視覺和圖像處理領域，用於目標檢測的特征描述器。這項技術是用來計算局部圖像梯度的方向信息的統計值。這種方法跟邊緣方向直方圖(edge orientation histograms)、尺度不變特征變換(scale-invariant feature transform descriptors) 以及形狀上下文方法( shape contexts)有很多相似之處，但與它們的不同點是：HOG描述器是在一個網格密集的大小統一的細胞單元(dense grid of uniformly spaced cells)上計算，而且為了提高性能，還采用了重疊的局部對比度歸一化(overlapping local contrast normalization)技術。
這篇文章的作者Navneet Dalal和Bill Triggs是法國國家計算機技術和控制研究所French National Institute for Research in?Computer Science and Control (INRIA)的研究員。他們在這篇文章中首次提出了HOG方法。這篇文章被發表在2005年的CVPR上。他們主要是將這種方法應用在靜態圖像中的行人 檢測上，但在后來，他們也將其應用在電影和視頻中的行人檢測，以及靜態圖像中的車輛和常見動物的檢測。
HOG描述器最重要的思想是：在一副 圖像中，局部目標的表象和形狀(appearance and shape)能夠被梯度或邊緣的方向密度分布很好地描述。具體的實現方法是：首先將圖像分成小的連通區域，我們把它叫細胞單元。然后采集細胞單元中各像素 點的梯度的或邊緣的方向直方圖。最后把這些直方圖組合起來就可以構成特征描述器。為了提高性能，我們還可以把這些局部直方圖在圖像的更大的范圍內(我們把 它叫區間或block)進行對比度歸一化(contrast-normalized)，所采用的方法是：先計算各直方圖在這個區間(block)中的密 度，然后根據這個密度對區間中的各個細胞單元做歸一化。通過這個歸一化后，能對光照變化和陰影獲得更好的效果。
與其他的特征描述方法相 比，HOG描述器后很多優點。首先，由於HOG方法是在圖像的局部細胞單元上操作，所以它對圖像幾何的(geometric)和光學的 (photometric)形變都能保持很好的不變性，這兩種形變只會出現在更大的空間領域上。其次，作者通過實驗發現，在粗的空域抽樣(coarse spatial sampling)、精細的方向抽樣(fine orientation sampling)以及較強的局部光學歸一化(strong local photometric normalization)等條件下，只要行人大體上能夠保持直立的姿勢，就容許行人有一些細微的肢體動作，這些細微的動作可以被忽略而不影響檢測效 果。綜上所述，HOG方法是特別適合於做圖像中的行人檢測的。
上圖是作者做的行人檢測試驗，其中(a)表示所有訓練圖像集 的平均梯度(average gradient across their training images)；(b)和(c)分別表示：圖像中每一個區間(block)上的最大最大正、負SVM權值；(d)表示一副測試圖像；(e)計算完R- HOG后的測試圖像；(f)和(g)分別表示被正、負SVM權值加權后的R-HOG圖像。
算法的實現：
色彩和伽馬歸一化 (color and gamma normalization)
作者分別在灰度空間、RGB色彩空間和LAB色彩空間上對圖像進行色彩和 伽馬歸一化，但實驗結果顯示，這個歸一化的預處理工作對最后的結果沒有影響，原因可能是：在后續步驟中也有歸一化的過程，那些過程可以取代這個預處理的歸 一化。所以，在實際應用中，這一步可以省略。
梯度的計算(Gradient computation)
最常用的方法是：簡單 地使用一個一維的離散微分模板(1-D centered point discrete derivative mask)在一個方向上或者同時在水平和垂直兩個方向上對圖像進行處理，更確切地說，這個方法需要使用下面的濾波器核濾除圖像中的色彩或變化劇烈的數據 (color or intensity data)
作者也嘗試了其他一些更復雜的模板，如3×3 Sobel 模板，或對角線模板(diagonal masks)，但是在這個行人檢測的實驗中，這些復雜模板的表現都較差，所以作者的結論是：模板越簡單，效果反而越好。作者也嘗試了在使用微分模板前加入 一個高斯平滑濾波，但是這個高斯平滑濾波的加入使得檢測效果更差，原因是：許多有用的圖像信息是來自變化劇烈的邊緣，而在計算梯度之前加入高斯濾波會把這 些邊緣濾除掉。
構建方向的直方圖(creating the orientation histograms)
第三步就是為 圖像的每個細胞單元構建梯度方向直方圖。細胞單元中的每一個像素點都為某個基於方向的直方圖通道(orientation-based histogram channel)投票。投票是采取加權投票(weighted voting)的方式，即每一票都是帶權值的，這個權值是根據該像素點的梯度幅度計算出來。可以采用幅值本身或者它的函數來表示這個權值，實際測試表明： 使用幅值來表示權值能獲得最佳的效果，當然，也可以選擇幅值的函數來表示，比如幅值的平方根(square root)、幅值的平方(square of the gradient magnitude)、幅值的截斷形式(clipped version of the magnitude)等。細胞單元可以是矩形的(rectangular)，也可以是星形的(radial)。直方圖通道是平均分布在0-1800(無 向)或0-3600(有向)范圍內。作者發現，采用無向的梯度和9個直方圖通道，能在行人檢測試驗中取得最佳的效果。
把細胞單元組 合成大的區間(grouping the cells together into larger blocks)
由於局部光照的變化 (variations of illumination)以及前景-背景對比度(foreground-background contrast)的變化，使得梯度強度(gradient strengths)的變化范圍非常大。這就需要對梯度強度做歸一化，作者采取的辦法是：把各個細胞單元組合成大的、空間上連通的區間(blocks)。 這樣以來，HOG描述器就變成了由各區間所有細胞單元的直方圖成分所組成的一個向量。這些區間是互有重疊的，這就意味著：每一個細胞單元的輸出都多次作用 於最終的描述器。區間有兩個主要的幾何形狀——矩形區間(R-HOG)和環形區間(C-HOG)。R-HOG區間大體上是一些方形的格子，它可以有三個參 數來表征：每個區間中細胞單元的數目、每個細胞單元中像素點的數目、每個細胞的直方圖通道數目。作者通過實驗表明，行人檢測的最佳參數設置是：3×3細胞 /區間、6×6像素/細胞、9個直方圖通道。作者還發現，在對直方圖做處理之前，給每個區間(block)加一個高斯空域窗口(Gaussian spatial window)是非常必要的，因為這樣可以降低邊緣的周圍像素點(pixels around the edge)的權重。
R- HOG跟SIFT描述器看起來很相似，但他們的不同之處是：R-HOG是在單一尺度下、密集的網格內、沒有對方向排序的情況下被計算出來(are computed in dense grids at some single scale without orientation alignment)；而SIFT描述器是在多尺度下、稀疏的圖像關鍵點上、對方向排序的情況下被計算出來(are computed at sparse scale-invariant key image points and are rotated to align orientation)。補充一點，R-HOG是各區間被組合起來用於對空域信息進行編碼(are used in conjunction to encode spatial form information)，而SIFT的各描述器是單獨使用的(are used singly)。
C- HOG區間(blocks)有兩種不同的形式，它們的區別在於：一個的中心細胞是完整的，一個的中心細胞是被分割的。如右圖所示：
作者發現 C-HOG的這兩種形式都能取得相同的效果。C-HOG區間(blocks)可以用四個參數來表征：角度盒子的個數(number of angular bins)、半徑盒子個數(number of radial bins)、中心盒子的半徑(radius of the center bin)、半徑的伸展因子(expansion factor for the radius)。通過實驗，對於行人檢測，最佳的參數設置為：4個角度盒子、2個半徑盒子、中心盒子半徑為4個像素、伸展因子為2。前面提到過，對於R- HOG，中間加一個高斯空域窗口是非常有必要的，但對於C-HOG，這顯得沒有必要。C-HOG看起來很像基於形狀上下文(Shape Contexts)的方法，但不同之處是：C-HOG的區間中包含的細胞單元有多個方向通道(orientation channels)，而基於形狀上下文的方法僅僅只用到了一個單一的邊緣存在數(edge presence count)。
區間歸一化 (Block normalization)
作者采用了四中不同的方法對區間進行歸一化，並對結果進行了比較。引入v表示一個還沒有被歸一 化的向量，它包含了給定區間(block)的所有直方圖信息。| | vk | |表示v的k階范數，這里的k去1、2。用e表示一個很小的常數。這時，歸一化因子可以表示如下：
L2-norm:
L1-norm:
L1-sqrt:
還 有第四種歸一化方式：L2-Hys，它可以通過先進行L2-norm，對結果進行截短(clipping)，然后再重新歸一化得到。作者發現：采用L2- Hys L2-norm 和 L1-sqrt方式所取得的效果是一樣的，L1-norm稍微表現出一點點不可靠性。但是對於沒有被歸一化的數據來說，這四種方法都表現出來顯著的改進。
SVM 分類器(SVM classifier)
最后一步就是把提取的HOG特征輸入到SVM分類器中，尋找一個最優超平面作為決策函數。作者采用 的方法是：使用免費的SVMLight軟件包加上HOG分類器來尋找測試圖像中的行人。
I=imread('05598.jpg');
c=imresize(I,[64?64],'bilinear');
f=hogcalculator(c,8,8,2,2,9,0.5,'localinterpolate','unsigned','l2hys');I=imread('05598.jpg');
c=imresize(I,[64 64],'bilinear');
f=hogcalculator(c,8,8,2,2,9,0.5,'localinterpolate','unsigned','l2hys');
函數舉例，將讀取圖像變成64*64，然后block=(8*2)*(8*2)=16*16 ，cell=8*8 ?bins=9
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?
%?set?default?parameters?value.
cellpw?=?2;
cellph?=?2;
nblockw?=?2;
nblockh?=?2;
nthet?=?9;
overlap?=?0.5;
isglobalinterpolate?=?'localinterpolate';
issigned?=?'unsigned';
normmethod?=?'l2hys';
else
if?nargin?
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
gradscalx?=?max(max(gradscalx(:,:,1),gradscalx(:,:,2)),?gradscalx(:,:,3));
gradscaly?=?max(max(gradscaly(:,:,1),gradscaly(:,:,2)),?gradscaly(:,:,3));
end
gradscal?=?sqrt(double(gradscalx.*gradscalx?+?gradscaly.*gradscaly));
%?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)?+?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?
gradorient(idx)?=?atan(gradscaly(idx)./gradscalxplus(idx))?+?pi;
idx?=?find(gradscalx?>=?0?&?gradscaly?
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(isglobalinterpolate,?'globalinterpolate')?==?1
ncellx?=?N?/?cellpw;
ncelly?=?M?/?cellph;
hist3dbig?=?zeros(ncelly+2,?ncellx+2,?nthet+2);
F?=?zeros(1,?(M/cellph-1)*(N/cellpw-1)*nblockw*nblockh*nthet);
glbalinter?=?1;
else
if?strcmp(isglobalinterpolate,?'localinterpolate')?==?1
hist3dbig?=?zeros(nblockh+2,?nblockw+2,?nthet+2);
F?=?zeros(1,?ntotalbh*ntotalbw*nblockw*nblockh*nthet);
glbalinter?=?0;
else
error('Incorrect?ISGLOBALINTERPOLATE?parameter.')
end
end
%?generate?the?matrix?for?storing?histogram?of?one?block;
sF?=?zeros(1,?nblockw*nblockh*nthet);
%?generate?gaussian?spatial?weight.
[gaussx,?gaussy]?=?meshgrid(0:(cellpw*nblockw-1),?0:(cellph*nblockh-1));
weight?=?exp(-((gaussx-(cellpw*nblockw-1)/2)...
.*(gaussx-(cellpw*nblockw-1)/2)+(gaussy-(cellph*nblockh-1)/2)...
.*(gaussy-(cellph*nblockh-1)/2))/(delta*delta));
%?vote?for?histogram.?there?are?two?situations?according?to?the?interpolate
%?condition('global'?interpolate?or?local?interpolate).?The?hist3d?which?is
%?generated?from?the?'bigger'?matrix?hist3dbig?is?the?final?histogram.
if?glbalinter?==?1
xbstep?=?nblockw*cellpw;
ybstep?=?nblockh*cellph;
else
xbstep?=?xbstride;
ybstep?=?ybstride;
end
%?block?slide?loop
for?btly?=?1:ybstep:ybstridend
for?btlx?=?1:xbstep:xbstridend
for?bi?=?1:(cellph*nblockh)
for?bj?=?1:(cellpw*nblockw)
i?=?btly?+?bi?-?1;
j?=?btlx?+?bj?-?1;
gaussweight?=?weight(bi,bj);
gs?=?gradscal(i,j);
go?=?gradorient(i,j);
if?glbalinter?==?1
jorbj?=?j;
iorbi?=?i;
else
jorbj?=?bj;
iorbi?=?bi;
end
%?calculate?bin?index?of?hist3dbig
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;
if?gs?==?0
continue;
end
binx2?=?binx1?+?1;
biny2?=?biny1?+?1;
binz2?=?binz1?+?1;
x1?=?(binx1-1.5)*cellpw?+?0.5;
y1?=?(biny1-1.5)*cellph?+?0.5;
z1?=?(binz1-1.5)*(or*pi/nthet);
%?trilinear?interpolation.
hist3dbig(biny1,binx1,binz1)?=...
hist3dbig(biny1,binx1,binz1)?+?gs*gaussweight...
*?(1-(jorbj-x1)/cellpw)*(1-(iorbi-y1)/cellph)...
*(1-(go-z1)/(or*pi/nthet));
hist3dbig(biny1,binx1,binz2)?=...
hist3dbig(biny1,binx1,binz2)?+?gs*gaussweight...
*?(1-(jorbj-x1)/cellpw)*(1-(iorbi-y1)/cellph)...
*((go-z1)/(or*pi/nthet));
hist3dbig(biny2,binx1,binz1)?=...
hist3dbig(biny2,binx1,binz1)?+?gs*gaussweight...
*?(1-(jorbj-x1)/cellpw)*((iorbi-y1)/cellph)...
*(1-(go-z1)/(or*pi/nthet));
hist3dbig(biny2,binx1,binz2)?=...
hist3dbig(biny2,binx1,binz2)?+?gs*gaussweight...
*?(1-(jorbj-x1)/cellpw)*((iorbi-y1)/cellph)...
*((go-z1)/(or*pi/nthet));
hist3dbig(biny1,binx2,binz1)?=...
hist3dbig(biny1,binx2,binz1)?+?gs*gaussweight...
*?((jorbj-x1)/cellpw)*(1-(iorbi-y1)/cellph)...
*(1-(go-z1)/(or*pi/nthet));
hist3dbig(biny1,binx2,binz2)?=...
hist3dbig(biny1,binx2,binz2)?+?gs*gaussweight...
*?((jorbj-x1)/cellpw)*(1-(iorbi-y1)/cellph)...
*((go-z1)/(or*pi/nthet));
hist3dbig(biny2,binx2,binz1)?=...
hist3dbig(biny2,binx2,binz1)?+?gs*gaussweight...
*?((jorbj-x1)/cellpw)*((iorbi-y1)/cellph)...
*(1-(go-z1)/(or*pi/nthet));
hist3dbig(biny2,binx2,binz2)?=...
hist3dbig(biny2,binx2,binz2)?+?gs*gaussweight...
*?((jorbj-x1)/cellpw)*((iorbi-y1)/cellph)...
*((go-z1)/(or*pi/nthet));
end
end
%?In?the?local?interpolate?condition.?F?is?generated?in?this?block
%?slide?loop.?hist3dbig?should?be?cleared?in?each?loop.
if?glbalinter?==?0
if?or?==?2
hist3dbig(:,:,2)?=?hist3dbig(:,:,2)...
+?hist3dbig(:,:,nthet+2);
hist3dbig(:,:,(nthet+1))?=...
hist3dbig(:,:,(nthet+1))?+?hist3dbig(:,:,1);
end
hist3d?=?hist3dbig(2:(nblockh+1),?2:(nblockw+1),?2:(nthet+1));
for?ibin?=?1:nblockh
for?jbin?=?1:nblockw
idsF?=?nthet*((ibin-1)*nblockw+jbin-1)+1;
idsF?=?idsF:(idsF+nthet-1);
sF(idsF)?=?hist3d(ibin,jbin,:);
end
end
iblock?=?((btly-1)/ybstride)*ntotalbw?+...
((btlx-1)/xbstride)?+?1;
idF?=?(iblock-1)*nblockw*nblockh*nthet+1;
idF?=?idF:(idF+nblockw*nblockh*nthet-1);
F(idF)?=?sF;
hist3dbig(:,:,:)?=?0;
end
end
end
%?In?the?global?interpolate?condition.?F?is?generated?here?outside?the
%?block?slide?loop
if?glbalinter?==?1
ncellx?=?N?/?cellpw;
ncelly?=?M?/?cellph;
if?or?==?2
hist3dbig(:,:,2)?=?hist3dbig(:,:,2)?+?hist3dbig(:,:,nthet+2);
hist3dbig(:,:,(nthet+1))?=?hist3dbig(:,:,(nthet+1))?+?hist3dbig(:,:,1);
end
hist3d?=?hist3dbig(2:(ncelly+1),?2:(ncellx+1),?2:(nthet+1));
iblock?=?1;
for?btly?=?1:ybstride:ybstridend
for?btlx?=?1:xbstride:xbstridend
binidx?=?floor((btlx-1)/cellpw)+1;
binidy?=?floor((btly-1)/cellph)+1;
idF?=?(iblock-1)*nblockw*nblockh*nthet+1;
idF?=?idF:(idF+nblockw*nblockh*nthet-1);
for?ibin?=?1:nblockh
for?jbin?=?1:nblockw
idsF?=?nthet*((ibin-1)*nblockw+jbin-1)+1;
idsF?=?idsF:(idsF+nthet-1);
sF(idsF)?=?hist3d(binidy+ibin-1,?binidx+jbin-1,?:);
end
end
F(idF)?=?sF;
iblock?=?iblock?+?1;
end
end
end
%?adjust?the?negative?value?caused?by?accuracy?of?floating-point
%?operations.these?value's?scale?is?very?small,?usually?at?E-03?magnitude
%?while?others?will?be?E+02?or?E+03?before?normalization.
F(F<0)?=?0;
%?block?normalization.
e?=?0.001;
l2hysthreshold?=?0.2;
fslidestep?=?nblockw*nblockh*nthet;
switch?normmethod
case?'none'
case?'l1'
for?fi?=?1:fslidestep:size(F,2)
div?=?sum(F(fi:(fi+fslidestep-1)));
F(fi:(fi+fslidestep-1))?=?F(fi:(fi+fslidestep-1))/(div+e);
end
case?'l1sqrt'
for?fi?=?1:fslidestep:size(F,2)
div?=?sum(F(fi:(fi+fslidestep-1)));
F(fi:(fi+fslidestep-1))?=?sqrt(F(fi:(fi+fslidestep-1))/(div+e));
end
case?'l2'
for?fi?=?1:fslidestep:size(F,2)
sF?=?F(fi:(fi+fslidestep-1)).*F(fi:(fi+fslidestep-1));
div?=?sqrt(sum(sF)+e*e);
F(fi:(fi+fslidestep-1))?=?F(fi:(fi+fslidestep-1))/div;
end
case?'l2hys'
for?fi?=?1:fslidestep:size(F,2)
sF?=?F(fi:(fi+fslidestep-1)).*F(fi:(fi+fslidestep-1));
div?=?sqrt(sum(sF)+e*e);
sF?=?F(fi:(fi+fslidestep-1))/div;
sF(sF>l2hysthreshold)?=?l2hysthreshold;
div?=?sqrt(sum(sF.*sF)+e*e);
F(fi:(fi+fslidestep-1))?=?sF/div;
end
otherwise
error('Incorrect?NORMMETHOD?parameter.');
end
%?the?following?code,?which?can?be?removed?because?of?having?no
%?contributions?to?HOG?feature?calculation,?are?just?for?result
%?demonstration?when?the?trilinear?interpolation?is?'global'?for?this
%?condition?is?easier?to?give?a?simple?and?intuitive?illustration.?so?in
%?'local'?condition?it?will?produce?error.
%?figure;
%?hold?on;
%?axis?equal;
%?xlim([0,?N]);
%?ylim([0,?M]);
%?for?u?=?1:(M/cellph)
%?????for?v?=?1:(N/cellpw)
%?????????cx?=?(v-1)*cellpw?+?cellpw/2?+?0.5;
%?????????cy?=?(u-1)*cellph?+?cellph/2?+?0.5;
%?????????hist3d(u,v,:)=0.9*min(cellpw,cellph)*hist3d(u,v,:)/max(hist3d(u,v,:));
%?????????for?t?=?1:nthet
%?????????????s?=?hist3d(u,v,t);
%?????????????thet?=?(t-1)*pi/nthet?+?pi*0.5/nthet;
%?????????????x1?=?cx?-?s*0.5*cos(thet);
%?????????????x2?=?cx?+?s*0.5*cos(thet);
%?????????????y1?=?cy?-?s*0.5*sin(thet);
%?????????????y2?=?cy?+?s*0.5*sin(thet);
%?????????????plot([x1,x2],[M-y1+1,M-y2+1]);
%?????????end
%?????end
%?endfunction 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 = 2;
cellph = 2;
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
gradscalx = max(max(gradscalx(:,:,1),gradscalx(:,:,2)), gradscalx(:,:,3));
gradscaly = max(max(gradscaly(:,:,1),gradscaly(:,:,2)), gradscaly(:,:,3));
end
gradscal = sqrt(double(gradscalx.*gradscalx + gradscaly.*gradscaly));
% 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) + 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.');
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(isglobalinterpolate, 'globalinterpolate') == 1
ncellx = N / cellpw;
ncelly = M / cellph;
hist3dbig = zeros(ncelly+2, ncellx+2, nthet+2);
F = zeros(1, (M/cellph-1)*(N/cellpw-1)*nblockw*nblockh*nthet);
glbalinter = 1;
else
if strcmp(isglobalinterpolate, 'localinterpolate') == 1
hist3dbig = zeros(nblockh+2, nblockw+2, nthet+2);
F = zeros(1, ntotalbh*ntotalbw*nblockw*nblockh*nthet);
glbalinter = 0;
else
error('Incorrect ISGLOBALINTERPOLATE parameter.')
end
end
% generate the matrix for storing histogram of one block;
sF = zeros(1, nblockw*nblockh*nthet);
% generate gaussian spatial weight.
[gaussx, gaussy] = meshgrid(0:(cellpw*nblockw-1), 0:(cellph*nblockh-1));
weight = exp(-((gaussx-(cellpw*nblockw-1)/2)...
.*(gaussx-(cellpw*nblockw-1)/2)+(gaussy-(cellph*nblockh-1)/2)...
.*(gaussy-(cellph*nblockh-1)/2))/(delta*delta));
% vote for histogram. there are two situations according to the interpolate
% condition('global' interpolate or local interpolate). The hist3d which is
% generated from the 'bigger' matrix hist3dbig is the final histogram.
if glbalinter == 1
xbstep = nblockw*cellpw;
ybstep = nblockh*cellph;
else
xbstep = xbstride;
ybstep = ybstride;
end
% block slide loop
for btly = 1:ybstep:ybstridend
for btlx = 1:xbstep:xbstridend
for bi = 1:(cellph*nblockh)
for bj = 1:(cellpw*nblockw)
i = btly + bi - 1;
j = btlx + bj - 1;
gaussweight = weight(bi,bj);
gs = gradscal(i,j);
go = gradorient(i,j);
if glbalinter == 1
jorbj = j;
iorbi = i;
else
jorbj = bj;
iorbi = bi;
end
% calculate bin index of hist3dbig
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;
if gs == 0
continue;
end
binx2 = binx1 + 1;
biny2 = biny1 + 1;
binz2 = binz1 + 1;
x1 = (binx1-1.5)*cellpw + 0.5;
y1 = (biny1-1.5)*cellph + 0.5;
z1 = (binz1-1.5)*(or*pi/nthet);
% trilinear interpolation.
hist3dbig(biny1,binx1,binz1) =...
hist3dbig(biny1,binx1,binz1) + gs*gaussweight...
* (1-(jorbj-x1)/cellpw)*(1-(iorbi-y1)/cellph)...
*(1-(go-z1)/(or*pi/nthet));
hist3dbig(biny1,binx1,binz2) =...
hist3dbig(biny1,binx1,binz2) + gs*gaussweight...
* (1-(jorbj-x1)/cellpw)*(1-(iorbi-y1)/cellph)...
*((go-z1)/(or*pi/nthet));
hist3dbig(biny2,binx1,binz1) =...
hist3dbig(biny2,binx1,binz1) + gs*gaussweight...
* (1-(jorbj-x1)/cellpw)*((iorbi-y1)/cellph)...
*(1-(go-z1)/(or*pi/nthet));
hist3dbig(biny2,binx1,binz2) =...
hist3dbig(biny2,binx1,binz2) + gs*gaussweight...
* (1-(jorbj-x1)/cellpw)*((iorbi-y1)/cellph)...
*((go-z1)/(or*pi/nthet));
hist3dbig(biny1,binx2,binz1) =...
hist3dbig(biny1,binx2,binz1) + gs*gaussweight...
* ((jorbj-x1)/cellpw)*(1-(iorbi-y1)/cellph)...
*(1-(go-z1)/(or*pi/nthet));
hist3dbig(biny1,binx2,binz2) =...
hist3dbig(biny1,binx2,binz2) + gs*gaussweight...
* ((jorbj-x1)/cellpw)*(1-(iorbi-y1)/cellph)...
*((go-z1)/(or*pi/nthet));
hist3dbig(biny2,binx2,binz1) =...
hist3dbig(biny2,binx2,binz1) + gs*gaussweight...
* ((jorbj-x1)/cellpw)*((iorbi-y1)/cellph)...
*(1-(go-z1)/(or*pi/nthet));
hist3dbig(biny2,binx2,binz2) =...
hist3dbig(biny2,binx2,binz2) + gs*gaussweight...
* ((jorbj-x1)/cellpw)*((iorbi-y1)/cellph)...
*((go-z1)/(or*pi/nthet));
end
end
% In the local interpolate condition. F is generated in this block
% slide loop. hist3dbig should be cleared in each loop.
if glbalinter == 0
if or == 2
hist3dbig(:,:,2) = hist3dbig(:,:,2)...
+ hist3dbig(:,:,nthet+2);
hist3dbig(:,:,(nthet+1)) =...
hist3dbig(:,:,(nthet+1)) + hist3dbig(:,:,1);
end
hist3d = hist3dbig(2:(nblockh+1), 2:(nblockw+1), 2:(nthet+1));
for ibin = 1:nblockh
for jbin = 1:nblockw
idsF = nthet*((ibin-1)*nblockw+jbin-1)+1;
idsF = idsF:(idsF+nthet-1);
sF(idsF) = hist3d(ibin,jbin,:);
end
end
iblock = ((btly-1)/ybstride)*ntotalbw +...
((btlx-1)/xbstride) + 1;
idF = (iblock-1)*nblockw*nblockh*nthet+1;
idF = idF:(idF+nblockw*nblockh*nthet-1);
F(idF) = sF;
hist3dbig(:,:,:) = 0;
end
end
end
% In the global interpolate condition. F is generated here outside the
% block slide loop
if glbalinter == 1
ncellx = N / cellpw;
ncelly = M / cellph;
if or == 2
hist3dbig(:,:,2) = hist3dbig(:,:,2) + hist3dbig(:,:,nthet+2);
hist3dbig(:,:,(nthet+1)) = hist3dbig(:,:,(nthet+1)) + hist3dbig(:,:,1);
end
hist3d = hist3dbig(2:(ncelly+1), 2:(ncellx+1), 2:(nthet+1));
iblock = 1;
for btly = 1:ybstride:ybstridend
for btlx = 1:xbstride:xbstridend
binidx = floor((btlx-1)/cellpw)+1;
binidy = floor((btly-1)/cellph)+1;
idF = (iblock-1)*nblockw*nblockh*nthet+1;
idF = idF:(idF+nblockw*nblockh*nthet-1);
for ibin = 1:nblockh
for jbin = 1:nblockw
idsF = nthet*((ibin-1)*nblockw+jbin-1)+1;
idsF = idsF:(idsF+nthet-1);
sF(idsF) = hist3d(binidy+ibin-1, binidx+jbin-1, :);
end
end
F(idF) = sF;
iblock = iblock + 1;
end
end
end
% adjust the negative value caused by accuracy of floating-point
% operations.these value's scale is very small, usually at E-03 magnitude
% while others will be E+02 or E+03 before normalization.
F(F<0) = 0;
% block normalization.
e = 0.001;
l2hysthreshold = 0.2;
fslidestep = nblockw*nblockh*nthet;
switch normmethod
case 'none'
case 'l1'
for fi = 1:fslidestep:size(F,2)
div = sum(F(fi:(fi+fslidestep-1)));
F(fi:(fi+fslidestep-1)) = F(fi:(fi+fslidestep-1))/(div+e);
end
case 'l1sqrt'
for fi = 1:fslidestep:size(F,2)
div = sum(F(fi:(fi+fslidestep-1)));
F(fi:(fi+fslidestep-1)) = sqrt(F(fi:(fi+fslidestep-1))/(div+e));
end
case 'l2'
for fi = 1:fslidestep:size(F,2)
sF = F(fi:(fi+fslidestep-1)).*F(fi:(fi+fslidestep-1));
div = sqrt(sum(sF)+e*e);
F(fi:(fi+fslidestep-1)) = F(fi:(fi+fslidestep-1))/div;
end
case 'l2hys'
for fi = 1:fslidestep:size(F,2)
sF = F(fi:(fi+fslidestep-1)).*F(fi:(fi+fslidestep-1));
div = sqrt(sum(sF)+e*e);
sF = F(fi:(fi+fslidestep-1))/div;
sF(sF>l2hysthreshold) = l2hysthreshold;
div = sqrt(sum(sF.*sF)+e*e);
F(fi:(fi+fslidestep-1)) = sF/div;
end
otherwise
error('Incorrect NORMMETHOD parameter.');
end
% the following code, which can be removed because of having no
% contributions to HOG feature calculation, are just for result
% demonstration when the trilinear interpolation is 'global' for this
% condition is easier to give a simple and intuitive illustration. so in
% 'local' condition it will produce error.
% figure;
% hold on;
% axis equal;
% xlim([0, N]);
% ylim([0, M]);
% for u = 1:(M/cellph)
% for v = 1:(N/cellpw)
% cx = (v-1)*cellpw + cellpw/2 + 0.5;
% cy = (u-1)*cellph + cellph/2 + 0.5;
% hist3d(u,v,:)=0.9*min(cellpw,cellph)*hist3d(u,v,:)/max(hist3d(u,v,:));
% for t = 1:nthet
% s = hist3d(u,v,t);
% thet = (t-1)*pi/nthet + pi*0.5/nthet;
% x1 = cx - s*0.5*cos(thet);
% x2 = cx + s*0.5*cos(thet);
% y1 = cy - s*0.5*sin(thet);
% y2 = cy + s*0.5*sin(thet);
% plot([x1,x2],[M-y1+1,M-y2+1]);
% end
% end
% end
總結
以上是生活随笔為你收集整理的matlab intergral,matlab學習：人臉識別之HOG（Histograms of Oriented Gradients）的全部內容，希望文章能夠幫你解決所遇到的問題。
- 上一篇： matlab中bp网络盲分离代码,利用m
- 下一篇： python mysql 编码方式,Py