hog訓練流程及原始碼分析
一、網上一些參考資料
在部落格目標檢測學習_1(用opencv自帶hog實現行人檢測) 中已經使用了opencv自帶的函式detectMultiScale()實現了對行人的檢測,當然了,該演算法採用的是hog演算法,那麼hog演算法是怎樣實現的呢?這一節就來簡單分析一下opencv中自帶 hog原始碼。
網上也有不少網友對opencv中的hog原始碼進行了分析,很不錯,看了很有收穫。比如:
http://blog.csdn.net/raocong2010/article/details/6239431
該部落格對該hog演算法中用到的block,cell等概念有一定的圖示解釋;
http://blog.csdn.net/pp5576155/article/details/7029699
該部落格是轉載的,裡面有opencv原始碼的一些註釋,很有幫助。
http://gz-ricky.blogbus.com/logs/85326280.html
本部落格對hog描述運算元長度的計算做了一定介紹。
http://hi.baidu.com/susongzhi/item/3a3c758d7ff5cbdc5e0ec172
該部落格對hog中快速演算法的三線插值將得很詳細。
http://blog.youtueye.com/work/opencv-hog-peopledetector-trainning.html
這篇部落格對hog怎樣訓練和檢測做了一定的講解。
二、關於原始碼的一些簡單說明
本文不是講解hog理論的,所以需要對hog演算法有一定了解,這些可以去參考hog提出者的博士論文,寫得很詳細。
按照正常流程,hog行人檢測分為訓練過程和檢測過程,訓練過程主要是訓練得到svm的係數。在opencv原始碼中直接採用訓練好了的svm係數,所以訓練過程原始碼中沒有涉及到多少。
首先還是對hog原始碼中一些固定引數來個簡單說明:
檢測視窗大小為128*64;
Block大小為16*16;
Cell大小為8*8;
Block在檢測視窗中上下移動尺寸為8*8;
1個cell的梯度直方圖化成9個bin;
滑動視窗在檢測圖片中滑動的尺寸為8*8;
程式碼中的一個hog描述子是針對一個檢測視窗而言的,所以一個檢測視窗共有105=((128-16)/8+1)*((64-16)/8+1)個block;一個block中有4個cell,而一個cell的hog描述子向量的長度為9;所以檢測視窗的hog向量長度=3780=105*4*9維。
三、hog訓練部分流程的簡單理解
雖然hog原始碼中很少涉及到訓練部分的程式碼,不過了解下訓練過程的流程會對整個檢測過程有個整體認識。
訓練過程中正樣本大小統一為128*64,即檢測視窗的大小;該樣本圖片可以包含1個或多個行人。對該圖片提前的hog特徵長度剛好為3780維,每一個特徵對應一個正樣本標籤進行訓練。在實際的訓練過程中,我們並不是去google上收集或者拍攝剛好128*64大小且有行人的圖片,而是收集包含行人的任意圖片(當然了,尺寸最好比128*64大),然後手工對這些正樣本進行標註,即對有行人的地方畫個矩形,其實也就是存了2個頂點的座標而已,並把這個矩形的資訊儲存起來;最好自己寫一個程式,每讀入一張圖片,就把矩形區域的內容截取出來並縮放到統一尺寸128*64,這樣,對處理過後的該圖片進行hog特徵提取就可以當做正樣本了。
負樣本不需要統一尺寸,只需比128*64大,且圖片中不能包含任何行人。實際過程中,由於是負樣本,裡面沒有目標資訊,所以不需要人工進行標註。程式中可以對該圖片隨機進行擷取128*64大小的圖片,並提取出其hog特徵作為負樣本。
四、hog行人檢測過程
檢測過程中採用的是滑動視窗法,對應本程式碼中,滑動視窗法的流程如下:
由上圖可以看出,檢測時,會對輸入圖片進行尺度縮放(一般是縮小),在每一層的影象上採用固定大小的滑動視窗(128*64)滑動,沒個滑動視窗都提取出hog特徵,送入到svm分類器中,看該視窗中是否有目標。有則存下目標區域來,無則繼續滑動。
檢測過程中用到的函式為detectMultiScale(),其引數分配圖如下:
五、計算檢測視窗中影象的梯度
計算梯度前如果需要gamma校正的話就先進行gamma校正,所謂的gamma校正就是把原來的每個通道畫素值範圍從0~255變換到0~15.97(255開根號)。據作者說這樣校正過後的影象計算的效果會更好,在計算梯度前不需要進行高斯濾波操作。
梯度的計算是分別計算水平梯度圖和垂直梯度圖,然後求幅值和相位。水平梯度卷積運算元為:
垂直梯度卷積運算元為:
在閱讀該原始碼的時候,要特別注意梯度幅值和角度的儲存方式。因為是對一個滑動窗口裡的影象進行的,所以梯度幅值和角度按照道理來說應該都是128*64=8192維的向量。但實際過程中這2者都是用的128*64*2=16384維的向量。為什麼呢?
因為這裡的梯度和角度都是用到了二線插值的。每一個點的梯度角度可能是0~180度之間的任意值,而程式中將其離散化為9個bin,即每個bin佔20度。所以滑動視窗中每個畫素點的梯度角度如果要離散化到這9個bin中,則一般它都會有2個相鄰的bin(如果恰好位於某個bin的中心,則可認為對該bin的權重為1即可)。從原始碼中可以看到梯度的幅值是用來計算梯度直方圖時權重投票的,所以每個畫素點的梯度幅值就分解到了其角度相鄰的2個bin了,越近的那個bin得到的權重越大。因此幅度影象用了2個通道,每個通道都是原畫素點幅度的一個分量。同理,不難理解,畫素點的梯度角度也用了2個通道,每個通道中儲存的是它相鄰2個bin的bin序號。序號小的放在第一通道。
二線插值的示意圖如下:
其中,假設那3條半徑為離散化後bin的中心,紅色虛線為畫素點O(畫素點在圓心處)的梯度方向,梯度幅值為A,該梯度方向與最近的相鄰bin為bin0,這兩者之間的夾角為a.這該畫素點O處儲存的梯度幅值第1通道為A*(1-a),第2通道為A*a;該畫素點O處儲存的角度第1通道為0(bin的序號為0),第2通道為1(bin的序號為1)。
另外在計算影象的梯度圖和相點陣圖時,如果該影象時3通道的,則3通道分別取梯度值,並且取梯度最大的那個通道的值為該點的梯度幅值。
六、HOG快取結構體
HOG快取思想是該程式作者加快hog演算法速度採用的一種記憶體優化技術。由於我們對每幅輸入圖片要進行4層掃描,分別為影象金字塔層,每層中滑動視窗,每個滑動視窗中滑動的block,每個block中的cell,其實還有每個cell中的畫素點;有這麼多層,每一層又是一個二維的,所以速度非常慢。作者的採用的思想是HOG快取,即把計算得到的每個滑動視窗的資料(其實最終是每個block的hog描述子向量)都存在記憶體查詢表中,由於滑動視窗在滑動時,很多個block都會重疊,因此重疊處計算過的block資訊就可以直接從查詢表中讀取,這樣就節省了很多時間。
在這個HOG儲存結構體中,會計算滑動視窗內的hog描述子,而這又涉及到滑動視窗,block,cell直接的關係,其之間的關係可以參考下面示意圖:
外面最大的為待檢測的圖片,對待檢測的圖片需要用滑動視窗進行滑動來判斷視窗中是否有目標,每個滑動視窗中又有很多個重疊移動的block,每個block中還有不重疊的cell。其實該程式的作者又將每個block中的畫素點對cell的貢獻不同,有將每個cell分成了4個區域,即圖中藍色虛線最小的框。
那麼block中不同的畫素點對它的cell(預設引數為1個block有4個cell)的影響是怎樣的呢?請看下面示意圖。
如果所示,黑色框代表1個block,紅實線隔開的為4個cell,每個cell用綠色虛線隔開的我們稱之為4個區域,所以該block中共有16個區域,分別為A、B、C、…、O、P。
程式中將這16個區域分為4組:
第1組:A、D、M、P;該組內的畫素點計算梯度方向直方圖時只對其所在的cell有貢獻。
第2組:B、C、N、O;該組內的畫素點計算梯度直方圖時對其所在的左右cell有貢獻。
第3組:E、I、H、L;該組內的畫素點計算梯度直方圖時對其所在的上下cell有貢獻。
第4組:F、G、J、K;該組內的畫素點對其上下左右的cell計算梯度直方圖時都有貢獻。
那到底是怎麼對cell貢獻的呢?舉個例子來說,E區域內的畫素點對cell0和cell2有貢獻。本來1個block對滑動視窗貢獻的向量維數為36維,即每個cell貢獻9維,其順序分別為cell0,cell1,cell2,cell3.而E區域內的畫素由於同時對cell0和cell2有貢獻,所以在計算E區域內的畫素梯度投票時,不僅要投向它本來的cell0,還要投向下面的cell2,即投向cell0和cell2有一個權重,該權重與該畫素點所在位置與cell0,cell2中心位置的距離有關。具體的關係可以去檢視原始碼。
該結構體變數記憶體分配圖如下,可以增強讀程式碼的直觀性:
在讀該部分原始碼時,需要特別注意以下幾個地方:
1) 結構體BlockData中有2個變數。1個BlockData結構體是對應的一個block資料。histOfs和imgOffset.其中histOfs表示為該block對整個滑動視窗內hog描述運算元的貢獻那部分向量的起始位置;imgOffset為該block在滑動視窗圖片中的座標(當然是指左上角座標)。
2) 結構體PixData中有5個變數,1個PixData結構體是對應的block中1個畫素點的資料。其中gradOfs表示該點的梯度幅度在滑動視窗圖片梯度幅度圖中的位置座標;qangleOfs表示該點的梯度角度在滑動視窗圖片梯度角度圖中的位置座標;histOfs[]表示該畫素點對1個或2個或4個cell貢獻的hog描述子向量的起始位置座標(比較抽象,需要看原始碼才懂)。histWeight[]表示該畫素點對1個或2個或4個cell貢獻的權重。gradWeight表示該點本身由於處在block中位置的不同因而對梯度直方圖貢獻也不同,其權值按照二維高斯分佈(以block中心為二維高斯的中心)來決定。
3) 程式中的count1,cout2,cout4分別表示該block中對1個cell、2個cell、4個cell有貢獻的畫素點的個數。
七、其他一些函式
該程式中還有一些其它的函式。
getblock()表示的是給定block在滑動視窗的位置以及圖片的hog快取指標,來獲得本次block中計算hog特徵所需要的資訊。
normalizeBlockHistogram()指對block獲取到的hog部分描述子進行歸一化,其實該歸一化有2層,具體看程式碼。
windowsInImage()實現的功能是給定測試圖片和滑動視窗移動的大小,來獲得該層中水平和垂直方向上需要滑動多少個滑動視窗。
getWindow()值獲得一個滑動視窗矩形。
compute()是實際上計算hog描述子的函式,在測試和訓練階段都能用到。
detect()是檢測目標是用到的函式,在detectMultiScale()函式內部被呼叫。
八、關於HOG的初始化
Hog初始化可以採用直接賦初值;也直接從檔案節點中讀取(有相應的格式,好像採用的是xml檔案格式);當然我們可以讀取初始值,也可以在程式中設定hog運算元的初始值並寫入檔案,這些工作可以採用原始碼中的read,write,load,save等函式來完成。
九、hog原始碼的註釋
在讀原始碼時,由於裡面用到了intel的ipp庫,優化了演算法的速度,所以在程式中遇到#ifdef HAVE_IPP後面的程式碼時,可以直接跳過不讀,直接讀#else後面的程式碼,這並不影響對原hog演算法的理解。
首先來看看hog原始碼中用到的標頭檔案目錄圖,如下:
下面是我對hog原始碼的一些註釋,由於本人接觸c++比較少,可能有些c++的語法常識也給註釋起來了,還望大家能理解。另外程式中還有一些細節沒有讀懂,或者說是註釋錯了的,大家可以一起來討論下,很多細節要在原始碼中才能看懂。
hog.cpp
/*M///////////////////////////////////////////////////////////////////////////////////////
//
// IMPORTANT: READ BEFORE DOWNLOADING, COPYING, INSTALLING OR USING.
//
// By downloading, copying, installing or using the software you agree to this license.
// If you do not agree to this license, do not download, install,
// copy or use the software.
//
//
// License Agreement
// For Open Source Computer Vision Library
//
// Copyright (C) 2000-2008, Intel Corporation, all rights reserved.
// Copyright (C) 2009, Willow Garage Inc., all rights reserved.
// Third party copyrights are property of their respective owners.
//
// Redistribution and use in source and binary forms, with or without modification,
// are permitted provided that the following conditions are met:
//
// * Redistribution's of source code must retain the above copyright notice,
// this list of conditions and the following disclaimer.
//
// * Redistribution's in binary form must reproduce the above copyright notice,
// this list of conditions and the following disclaimer in the documentation
// and/or other materials provided with the distribution.
//
// * The name of the copyright holders may not be used to endorse or promote products
// derived from this software without specific prior written permission.
//
// This software is provided by the copyright holders and contributors "as is" and
// any express or implied warranties, including, but not limited to, the implied
// warranties of merchantability and fitness for a particular purpose are disclaimed.
// In no event shall the Intel Corporation or contributors be liable for any direct,
// indirect, incidental, special, exemplary, or consequential damages
// (including, but not limited to, procurement of substitute goods or services;
// loss of use, data, or profits; or business interruption) however caused
// and on any theory of liability, whether in contract, strict liability,
// or tort (including negligence or otherwise) arising in any way out of
// the use of this software, even if advised of the possibility of such damage.
//
//M*/
#include "precomp.hpp"
#include <iterator>
#ifdef HAVE_IPP
#include "ipp.h"
#endif
/****************************************************************************************\
The code below is implementation of HOG (Histogram-of-Oriented Gradients)
descriptor and object detection, introduced by Navneet Dalal and Bill Triggs.
The computed feature vectors are compatible with the
INRIA Object Detection and Localization Toolkit
(http://pascal.inrialpes.fr/soft/olt/)
\****************************************************************************************/
namespace cv
{
size_t HOGDescriptor::getDescriptorSize() const
{
//下面2個語句是保證block中有整數個cell;保證block在視窗中能移動整數次
CV_Assert(blockSize.width % cellSize.width == 0 &&
blockSize.height % cellSize.height == 0);
CV_Assert((winSize.width - blockSize.width) % blockStride.width == 0 &&
(winSize.height - blockSize.height) % blockStride.height == 0 );
//返回的nbins是每個視窗中檢測到的hog向量的維數
return (size_t)nbins*
(blockSize.width/cellSize.width)*
(blockSize.height/cellSize.height)*
((winSize.width - blockSize.width)/blockStride.width + 1)*
((winSize.height - blockSize.height)/blockStride.height + 1);
}
//winSigma到底是什麼作用呢?
double HOGDescriptor::getWinSigma() const
{
return winSigma >= 0 ? winSigma : (blockSize.width + blockSize.height)/8.;
}
//svmDetector是HOGDescriptor內的一個成員變數,資料型別為向量vector。
//用來儲存hog特徵用於svm分類時的係數的.
//該函式返回為真的實際含義是什麼呢?保證與hog特徵長度相同,或者相差1,但為什麼
//相差1也可以呢?
bool HOGDescriptor::checkDetectorSize() const
{
size_t detectorSize = svmDetector.size(), descriptorSize = getDescriptorSize();
return detectorSize == 0 ||
detectorSize == descriptorSize ||
detectorSize == descriptorSize + 1;
}
void HOGDescriptor::setSVMDetector(InputArray _svmDetector)
{
//這裡的convertTo函式只是將影象Mat屬性更改,比如說通道數,矩陣深度等。
//這裡是將輸入的svm係數矩陣全部轉換成浮點型。
_svmDetector.getMat().convertTo(svmDetector, CV_32F);
CV_Assert( checkDetectorSize() );
}
#define CV_TYPE_NAME_HOG_DESCRIPTOR "opencv-object-detector-hog"
//FileNode是opencv的core中的一個檔案儲存節點類,這個節點用來儲存讀取到的每一個檔案元素。
//一般是讀取XML和YAML格式的檔案
//又因為該函式是把檔案節點中的內容讀取到其類的成員變數中,所以函式後面不能有關鍵字const
bool HOGDescriptor::read(FileNode& obj)
{
//isMap()是用來判斷這個節點是不是一個對映型別,如果是對映型別,則每個節點都與
//一個名字對應起來。因此這裡的if語句的作用就是需讀取的檔案node是一個對映型別
if( !obj.isMap() )
return false;
//中括號中的"winSize"是指返回名為winSize的一個節點,因為已經知道這些節點是mapping型別
//也就是說都有一個對應的名字。
FileNodeIterator it = obj["winSize"].begin();
//操作符>>為從節點中讀入資料,這裡是將it指向的節點資料依次讀入winSize.width,winSize.height
//下面的幾條語句功能類似
it >> winSize.width >> winSize.height;
it = obj["blockSize"].begin();
it >> blockSize.width >> blockSize.height;
it = obj["blockStride"].begin();
it >> blockStride.width >> blockStride.height;
it = obj["cellSize"].begin();
it >> cellSize.width >> cellSize.height;
obj["nbins"] >> nbins;
obj["derivAperture"] >> derivAperture;
obj["winSigma"] >> winSigma;
obj["histogramNormType"] >> histogramNormType;
obj["L2HysThreshold"] >> L2HysThreshold;
obj["gammaCorrection"] >> gammaCorrection;
obj["nlevels"] >> nlevels;
//isSeq()是判斷該節點內容是不是一個序列
FileNode vecNode = obj["SVMDetector"];
if( vecNode.isSeq() )
{
vecNode >> svmDetector;
CV_Assert(checkDetectorSize());
}
//上面的都讀取完了後就返回讀取成功標誌
return true;
}
void HOGDescriptor::write(FileStorage& fs, const String& objName) const
{
//將objName名字輸入到檔案fs中
if( !objName.empty() )
fs << objName;
fs << "{" CV_TYPE_NAME_HOG_DESCRIPTOR
//下面幾句依次將hog描述子內的變數輸入到檔案fs中,且每次輸入前都輸入
//一個名字與其對應,因此這些節點是mapping型別。
<< "winSize" << winSize
<< "blockSize" << blockSize
<< "blockStride" << blockStride
<< "cellSize" << cellSize
<< "nbins" << nbins
<< "derivAperture" << derivAperture
<< "winSigma" << getWinSigma()
<< "histogramNormType" << histogramNormType
<< "L2HysThreshold" << L2HysThreshold
<< "gammaCorrection" << gammaCorrection
<< "nlevels" << nlevels;
if( !svmDetector.empty() )
//svmDetector則是直接輸入序列,也有對應的名字。
fs << "SVMDetector" << "[:" << svmDetector << "]";
fs << "}";
}
//從給定的檔案中讀取引數
bool HOGDescriptor::load(const String& filename, const String& objname)
{
FileStorage fs(filename, FileStorage::READ);
//一個檔案節點有很多葉子,所以一個檔案節點包含了很多內容,這裡當然是包含的
//HOGDescriptor需要的各種引數了。
FileNode obj = !objname.empty() ? fs[objname] : fs.getFirstTopLevelNode();
return read(obj);
}
//將類中的引數以檔案節點的形式寫入檔案中。
void HOGDescriptor::save(const String& filename, const String& objName) const
{
FileStorage fs(filename, FileStorage::WRITE);
write(fs, !objName.empty() ? objName : FileStorage::getDefaultObjectName(filename));
}
//複製HOG描述子到c中
void HOGDescriptor::copyTo(HOGDescriptor& c) const
{
c.winSize = winSize;
c.blockSize = blockSize;
c.blockStride = blockStride;
c.cellSize = cellSize;
c.nbins = nbins;
c.derivAperture = derivAperture;
c.winSigma = winSigma;
c.histogramNormType = histogramNormType;
c.L2HysThreshold = L2HysThreshold;
c.gammaCorrection = gammaCorrection;
//vector型別也可以用等號賦值
c.svmDetector = svmDetector; c.nlevels = nlevels; }
//計算影象img的梯度幅度影象grad和梯度方向影象qangle.
//paddingTL為需要在原影象img左上角擴增的尺寸,同理paddingBR
//為需要在img影象右下角擴增的尺寸。
void HOGDescriptor::computeGradient(const Mat& img, Mat& grad, Mat& qangle,
Size paddingTL, Size paddingBR) const
{
//該函式只能計算8位整型深度的單通道或者3通道影象.
CV_Assert( img.type() == CV_8U || img.type() == CV_8UC3 );
//將影象按照輸入引數進行擴充,這裡不是為了計算邊緣梯度而做的擴充,因為
//為了邊緣梯度而擴充是在後面的程式碼完成的,所以這裡為什麼擴充暫時還不明白。
Size gradsize(img.cols + paddingTL.width + paddingBR.width,
img.rows + paddingTL.height + paddingBR.height);
grad.create(gradsize, CV_32FC2); // <magnitude*(1-alpha), magnitude*alpha>
qangle.create(gradsize, CV_8UC2); // [0..nbins-1] - quantized gradient orientation
Size wholeSize;
Point roiofs;
//locateROI在此處是如果img影象是從其它父影象中某一部分得來的,那麼其父影象
//的大小尺寸就為wholeSize了,img影象左上角相對於父影象的位置點就為roiofs了。
//對於正樣本,其父影象就是img了,所以這裡的wholeSize就和img.size()是一樣的,
//對應負樣本,這2者不同;因為裡面的關係比較不好懂,這裡權且將wholesSize理解為
//img的size,所以roiofs就應當理解為Point(0, 0)了。
img.locateROI(wholeSize, roiofs);
int i, x, y;
int cn = img.channels();
//_lut為行向量,用來作為浮點畫素值的儲存查詢表
Mat_<float> _lut(1, 256);
const float* lut = &_lut(0,0);
//gamma校正指的是將0~256的畫素值全部開根號,即範圍縮小了,且變換範圍都不成線性了,
if( gammaCorrection )
for( i = 0; i < 256; i++ )
_lut(0,i) = std::sqrt((float)i);
else
for( i = 0; i < 256; i++ )
_lut(0,i) = (float)i;
//建立長度為gradsize.width+gradsize.height+4的整型buffer
AutoBuffer<int> mapbuf(gradsize.width + gradsize.height + 4);
int* xmap = (int*)mapbuf + 1;
int* ymap = xmap + gradsize.width + 2;
//言外之意思borderType就等於4了,因為opencv的原始碼中是如下定義的。
//#define IPL_BORDER_REFLECT_101 4
//enum{...,BORDER_REFLECT_101=IPL_BORDER_REFLECT_101,...}
//borderType為邊界擴充後所填充畫素點的方式。
/*
Various border types, image boundaries are denoted with '|'
* BORDER_REPLICATE: aaaaaa|abcdefgh|hhhhhhh
* BORDER_REFLECT: fedcba|abcdefgh|hgfedcb
* BORDER_REFLECT_101: gfedcb|abcdefgh|gfedcba
* BORDER_WRAP: cdefgh|abcdefgh|abcdefg
* BORDER_CONSTANT: iiiiii|abcdefgh|iiiiiii with some specified 'i'
*/
const int borderType = (int)BORDER_REFLECT_101;
for( x = -1; x < gradsize.width + 1; x++ )
/*int borderInterpolate(int p, int len, int borderType)
其中引數p表示的是擴充後圖像的一個座標,相對於對應的座標軸而言;
len引數表示對應源影象的一個座標軸的長度;borderType為擴充型別,
在上面已經有過介紹.
所以這個函式的作用是從擴充後的畫素點座標推斷出源影象中對應該點
的座標值。
*/
//這裡的xmap和ymap實際含義是什麼呢?其實xmap向量裡面存的就是
//擴充後圖像第一行畫素點對應與原影象img中的畫素橫座標,可以看
//出,xmap向量中有些元素的值是相同的,因為擴充影象肯定會對應
//到原影象img中的某一位置,而img本身尺寸內的畫素也會對應該位置。
//同理,ymap向量裡面存的是擴充後圖像第一列畫素點對應於原圖想img
//中的畫素縱座標。
xmap[x] = borderInterpolate(x - paddingTL.width + roiofs.x,
wholeSize.width, borderType) - roiofs.x;
for( y = -1; y < gradsize.height + 1; y++ )
ymap[y] = borderInterpolate(y - paddingTL.height + roiofs.y,
wholeSize.height, borderType) - roiofs.y;
// x- & y- derivatives for the whole row
int width = gradsize.width;
AutoBuffer<float> _dbuf(width*4);
float* dbuf = _dbuf;
//DX為水平梯度圖,DY為垂直梯度圖,Mag為梯度幅度圖,Angle為梯度角度圖
//該構造方法的第4個引數表示矩陣Mat的資料在記憶體中存放的位置。由此可以
//看出,這4幅影象在記憶體中是連續儲存的。
Mat Dx(1, width, CV_32F, dbuf);
Mat Dy(1, width, CV_32F, dbuf + width);
Mat Mag(1, width, CV_32F, dbuf + width*2);
Mat Angle(1, width, CV_32F, dbuf + width*3);
int _nbins = nbins;
//angleScale==9/pi;
float angleScale = (float)(_nbins/CV_PI);
#ifdef HAVE_IPP
Mat lutimg(img.rows,img.cols,CV_MAKETYPE(CV_32F,cn));
Mat hidxs(1, width, CV_32F);
Ipp32f* pHidxs = (Ipp32f*)hidxs.data;
Ipp32f* pAngles = (Ipp32f*)Angle.data;
IppiSize roiSize;
roiSize.width = img.cols;
roiSize.height = img.rows;
for( y = 0; y < roiSize.height; y++ )
{
const uchar* imgPtr = img.data + y*img.step;
float* imglutPtr = (float*)(lutimg.data + y*lutimg.step);
for( x = 0; x < roiSize.width*cn; x++ )
{
imglutPtr[x] = lut[imgPtr[x]];
}
}
#endif
for( y = 0; y < gradsize.height; y++ )
{
#ifdef HAVE_IPP
const float* imgPtr = (float*)(lutimg.data + lutimg.step*ymap[y]);
const float* prevPtr = (float*)(lutimg.data + lutimg.step*ymap[y-1]);
const float* nextPtr = (float*)(lutimg.data + lutimg.step*ymap[y+1]);
#else
//imgPtr在這裡指的是img影象的第y行首地址;prePtr指的是img第y-1行首地址;
//nextPtr指的是img第y+1行首地址;
const uchar* imgPtr = img.data + img.step*ymap[y];
const uchar* prevPtr = img.data + img.step*ymap[y-1];
const uchar* nextPtr = img.data + img.step*ymap[y+1];
#endif
float* gradPtr = (float*)grad.ptr(y);
uchar* qanglePtr = (uchar*)qangle.ptr(y);
//輸入影象img為單通道影象時的計算
if( cn == 1 )
{
for( x = 0; x < width; x++ )
{
int x1 = xmap[x];
#ifdef HAVE_IPP
dbuf[x] = (float)(imgPtr[xmap[x+1]] - imgPtr[xmap[x-1]]);
dbuf[width + x] = (float)(nextPtr[x1] - prevPtr[x1]);
#else
//下面2句把Dx,Dy就計算出來了,因為其對應的記憶體都在dbuf中
dbuf[x] = (float)(lut[imgPtr[xmap[x+1]]] - lut[imgPtr[xmap[x-1]]]);
dbuf[width + x] = (float)(lut[nextPtr[x1]] - lut[prevPtr[x1]]);
#endif
}
}
//當cn==3時,也就是輸入影象為3通道影象時的處理。
else
{
for( x = 0; x < width; x++ )
{
//x1表示第y行第x1列的地址
int x1 = xmap[x]*3;
float dx0, dy0, dx, dy, mag0, mag;
#ifdef HAVE_IPP
const float* p2 = imgPtr + xmap[x+1]*3;
const float* p0 = imgPtr + xmap[x-1]*3;
dx0 = p2[2] - p0[2];
dy0 = nextPtr[x1+2] - prevPtr[x1+2];
mag0 = dx0*dx0 + dy0*dy0;
dx = p2[1] - p0[1];
dy = nextPtr[x1+1] - prevPtr[x1+1];
mag = dx*dx + dy*dy;
if( mag0 < mag )
{
dx0 = dx;
dy0 = dy;
mag0 = mag;
}
dx = p2[0] - p0[0];
dy = nextPtr[x1] - prevPtr[x1];
mag = dx*dx + dy*dy;
#else
//p2為第y行第x+1列的地址
//p0為第y行第x-1列的地址
const uchar* p2 = imgPtr + xmap[x+1]*3;
const uchar* p0 = imgPtr + xmap[x-1]*3;
//計算第2通道的幅值
dx0 = lut[p2[2]] - lut[p0[2]];
dy0 = lut[nextPtr[x1+2]] - lut[prevPtr[x1+2]];
mag0 = dx0*dx0 + dy0*dy0;
//計算第1通道的幅值
dx = lut[p2[1]] - lut[p0[1]];
dy = lut[nextPtr[x1+1]] - lut[prevPtr[x1+1]];
mag = dx*dx + dy*dy;
//取幅值最大的那個通道
if( mag0 < mag )
{
dx0 = dx;
dy0 = dy;
mag0 = mag;
}
//計算第0通道的幅值
dx = lut[p2[0]] - lut[p0[0]];
dy = lut[nextPtr[x1]] - lut[prevPtr[x1]];
mag = dx*dx + dy*dy;
#endif
//取幅值最大的那個通道
if( mag0 < mag )
{
dx0 = dx;
dy0 = dy;
mag0 = mag;
}
//最後求出水平和垂直方向上的梯度影象
dbuf[x] = dx0;
dbuf[x+width] = dy0;
}
}
#ifdef HAVE_IPP
ippsCartToPolar_32f((const Ipp32f*)Dx.data, (const Ipp32f*)Dy.data, (Ipp32f*)Mag.data, pAngles, width);
for( x = 0; x < width; x++ )
{
if(pAngles[x] < 0.f)
pAngles[x] += (Ipp32f)(CV_PI*2.);
}
ippsNormalize_32f(pAngles, pAngles, width, 0.5f/angleScale, 1.f/angleScale);
ippsFloor_32f(pAngles,(Ipp32f*)hidxs.data,width);
ippsSub_32f_I((Ipp32f*)hidxs.data,pAngles,width);
ippsMul_32f_I((Ipp32f*)Mag.data,pAngles,width);
ippsSub_32f_I(pAngles,(Ipp32f*)Mag.data,width);
ippsRealToCplx_32f((Ipp32f*)Mag.data,pAngles,(Ipp32fc*)gradPtr,width);
#else
//cartToPolar()函式是計算2個矩陣對應元素的幅度和角度,最後一個引數為是否
//角度使用度數表示,這裡為false表示不用度數表示,即用弧度表示。
//如果只需計算2個矩陣對應元素的幅度影象,可以採用magnitude()函式。
//-pi/2<Angle<pi/2;
cartToPolar( Dx, Dy, Mag, Angle, false );
#endif
for( x = 0; x < width; x++ )
{
#ifdef HAVE_IPP
int hidx = (int)pHidxs[x];
#else
//-5<angle<4
float mag = dbuf[x+width*2], angle = dbuf[x+width*3]*angleScale - 0.5f;
//cvFloor()返回不大於引數的最大整數
//hidx={-5,-4,-3,-2,-1,0,1,2,3,4};
int hidx = cvFloor(angle);
//0<=angle<1;angle表示的意思是與其相鄰的較小的那個bin的弧度距離(即弧度差)
angle -= hidx;
//gradPtr為grad影象的指標
//gradPtr[x*2]表示的是與x處梯度方向相鄰較小的那個bin的幅度權重;
//gradPtr[x*2+1]表示的是與x處梯度方向相鄰較大的那個bin的幅度權重
gradPtr[x*2] = mag*(1.f - angle);
gradPtr[x*2+1] = mag*angle;
#endif
if( hidx < 0 )
hidx += _nbins;
else if( hidx >= _nbins )
hidx -= _nbins;
assert( (unsigned)hidx < (unsigned)_nbins );
qanglePtr[x*2] = (uchar)hidx;
hidx++;
//-1在補碼中的表示為11111111,與-1相與的話就是自己本身了;
//0在補碼中的表示為00000000,與0相與的結果就是0了.
hidx &= hidx < _nbins ? -1 : 0;
qanglePtr[x*2+1] = (uchar)hidx;
}
}
}
struct HOGCache
{
struct BlockData
{
BlockData() : histOfs(0), imgOffset() {}
int histOfs;
Point imgOffset;
};
struct PixData
{
size_t gradOfs, qangleOfs;
int histOfs[4];
float histWeights[4];
float gradWeight;
};
HOGCache();
HOGCache(const HOGDescriptor* descriptor,
const Mat& img, Size paddingTL, Size paddingBR,
bool useCache, Size cacheStride);
virtual ~HOGCache() {};
virtual void init(const HOGDescriptor* descriptor,
const Mat& img, Size paddingTL, Size paddingBR,
bool useCache, Size cacheStride);
Size windowsInImage(Size imageSize, Size winStride) const;
Rect getWindow(Size imageSize, Size winStride, int idx) const;
const float* getBlock(Point pt, float* buf);
virtual void normalizeBlockHistogram(float* histogram) const;
vector<PixData> pixData;
vector<BlockData> blockData;
bool useCache;
vector<int> ymaxCached;
Size winSize, cacheStride;
Size nblocks, ncells;
int blockHistogramSize;
int count1, count2, count4;
Point imgoffset;
Mat_<float> blockCache;
Mat_<uchar> blockCacheFlags;
Mat grad, qangle;
const HOGDescriptor* descriptor;
};
//預設的建構函式,不使用cache,塊的直方圖向量大小為0等
HOGCache::HOGCache()
{
useCache = false;
blockHistogramSize = count1 = count2 = count4 = 0;
descriptor = 0;
}
//帶參的初始化函式,採用內部的init函式進行初始化
HOGCache::HOGCache(const HOGDescriptor* _descriptor,
const Mat& _img, Size _paddingTL, Size _paddingBR,
bool _useCache, Size _cacheStride)
{
init(_descriptor, _img, _paddingTL, _paddingBR, _useCache, _cacheStride);
}
//HOGCache結構體的初始化函式
void HOGCache::init(const HOGDescriptor* _descriptor,
const Mat& _img, Size _paddingTL, Size _paddingBR,
bool _useCache, Size _cacheStride)
{
descriptor = _descriptor;
cacheStride = _cacheStride;
useCache = _useCache;
//首先呼叫computeGradient()函式計算輸入影象的權值梯度幅度圖和角度量化圖
descriptor->computeGradient(_img, grad, qangle, _paddingTL, _paddingBR);
//imgoffset是Point型別,而_paddingTL是Size型別,雖然型別不同,但是2者都是
//一個二維座標,所以是在opencv中是允許直接賦值的。
imgoffset = _paddingTL;
winSize = descriptor->winSize;
Size blockSize = descriptor->blockSize;
Size blockStride = descriptor->blockStride;
Size cellSize = descriptor->cellSize;
int i, j, nbins = descriptor->nbins;
//rawBlockSize為block中包含畫素點的個數
int rawBlockSize = blockSize.width*blockSize.height;
//nblocks為Size型別,其長和寬分別表示一個視窗中水平方向和垂直方向上block的
//個數(需要考慮block在視窗中的移動)
nblocks = Size((winSize.width - blockSize.width)/blockStride.width + 1,
(winSize.height - blockSize.height)/blockStride.height + 1);
//ncells也是Size型別,其長和寬分別表示一個block中水平方向和垂直方向容納下
//的cell個數
ncells = Size(blockSize.width/cellSize.width, blockSize.height/cellSize.height);
//blockHistogramSize表示一個block中貢獻給hog描述子向量的長度
blockHistogramSize = ncells.width*ncells.height*nbins;
if( useCache )
{
//cacheStride= _cacheStride,即其大小是由引數傳入的,表示的是視窗移動的大小
//cacheSize長和寬表示擴充後的影象cache中,block在水平方向和垂直方向出現的個數
Size cacheSize((grad.cols - blockSize.width)/cacheStride.width+1,
(winSize.height/cacheStride.height)+1);
//blockCache為一個float型的Mat,注意其列數的值
blockCache.create(cacheSize.height, cacheSize.width*blockHistogramSize);
//blockCacheFlags為一個uchar型的Mat
blockCacheFlags.create(cacheSize);
size_t cacheRows = blockCache.rows;
//ymaxCached為vector<int>型別
//Mat::resize()為矩陣的一個方法,只是改變矩陣的行數,與單獨的resize()函式不相同。
ymaxCached.resize(cacheRows);
//ymaxCached向量內部全部初始化為-1
for(size_t ii = 0; ii < cacheRows; ii++ )
ymaxCached[ii] = -1;
}
//weights為一個尺寸為blockSize的二維高斯表,下面的程式碼就是計算二維高斯的係數
Mat_<float> weights(blockSize);
float sigma = (float)descriptor->getWinSigma();
float scale = 1.f/(sigma*sigma*2);
for(i = 0; i < blockSize.height; i++)
for(j = 0; j < blockSize.width; j++)
{
float di = i - blockSize.height*0.5f;
float dj = j - blockSize.width*0.5f;
weights(i,j) = std::exp(-(di*di + dj*dj)*scale);
}
//vector<BlockData> blockData;而BlockData為HOGCache的一個結構體成員
//nblocks.width*nblocks.height表示一個檢測視窗中block的個數,
//而cacheSize.width*cacheSize.heigh表示一個已經擴充的圖片中的block的個數
blockData.resize(nblocks.width*nblocks.height);
//vector<PixData> pixData;同理,Pixdata也為HOGCache中的一個結構體成員
//rawBlockSize表示每個block中畫素點的個數
//resize表示將其轉換成列向量
pixData.resize(rawBlockSize*3);
// Initialize 2 lookup tables, pixData & blockData.
// Here is why:
//
// The detection algorithm runs in 4 nested loops (at each pyramid layer):
// loop over the windows within the input image
// loop over the blocks within each window
// loop over the cells within each block
// loop over the pixels in each cell
//
// As each of the loops runs over a 2-dimensional array,
// we could get 8(!) nested loops in total, which is very-very slow.
//
// To speed the things up, we do the following:
// 1. loop over windows is unrolled in the HOGDescriptor::{compute|detect} methods;
// inside we compute the current search window using getWindow() method.
// Yes, it involves some overhead (function call + couple of divisions),
// but it's tiny in fact.
// 2. loop over the blocks is also unrolled. Inside we use pre-computed blockData[j]
// to set up gradient and histogram pointers.
// 3. loops over cells and pixels in each cell are merged
// (since there is no overlap between cells, each pixel in the block is processed once)
// and also unrolled. Inside we use PixData[k] to access the gradient values and
// update the histogram
//count1,count2,count4分別表示block中同時對1個cell,2個cell,4個cell有貢獻的畫素點的個數。
count1 = count2 = count4 = 0;
for( j = 0; j < blockSize.width; j++ )
for( i = 0; i < blockSize.height; i++ )
{
PixData* data = 0;
//cellX和cellY表示的是block內該畫素點所在的cell橫座標和縱座標索引,以小數的形式存在。
float cellX = (j+0.5f)/cellSize.width - 0.5f;
float cellY = (i+0.5f)/cellSize.height - 0.5f;
//cvRound返回最接近引數的整數;cvFloor返回不大於引數的整數;cvCeil返回不小於引數的整數
//icellX0和icellY0表示所在cell座標索引,索引值為該畫素點相鄰cell的那個較小的cell索引
//當然此處就是由整數的形式存在了。
//按照預設的係數的話,icellX0和icellY0只可能取值-1,0,1,且當i和j<3.5時對應的值才取-1
//當i和j>11.5時取值為1,其它時刻取值為0(注意i,j最大是15,從0開始的)
int icellX0 = cvFloor(cellX);
int icellY0 = cvFloor(cellY);
int icellX1 = icellX0 + 1, icellY1 = icellY0 + 1;
//此處的cellx和celly表示的是真實索引值與最近鄰cell索引值之間的差,
//為後面計算同一畫素對不同cell中的hist權重的計算。
cellX -= icellX0;
cellY -= icellY0;
//滿足這個if條件說明icellX0只能為0,也就是說block橫座標在(3.5,11.5)之間時
if( (unsigned)icellX0 < (unsigned)ncells.width &&
(unsigned)icellX1 < (unsigned)ncells.width )
{
//滿足這個if條件說明icellY0只能為0,也就是說block縱座標在(3.5,11.5)之間時
if( (unsigned)icellY0 < (unsigned)ncells.height &&
(unsigned)icellY1 < (unsigned)ncells.height )
{
//同時滿足上面2個if語句的畫素對4個cell都有權值貢獻
//rawBlockSize表示的是1個block中儲存畫素點的個數
//而pixData的尺寸大小為block中畫素點的3倍,其定義如下:
//pixData.resize(rawBlockSize*3);
//pixData的前面block畫素大小的記憶體為儲存只對block中一個cell
//有貢獻的pixel;中間block畫素大小的記憶體儲存對block中同時2個
//cell有貢獻的pixel;最後面的為對block中同時4個cell都有貢獻
//的pixel
data = &pixData[rawBlockSize*2 + (count4++)];
//下面計算出的結果為0
data->histOfs[0] = (icellX0*ncells.height + icellY0)*nbins;
//為該畫素點對cell0的權重
data->histWeights[0] = (1.f - cellX)*(1.f - cellY);
//下面計算出的結果為18
data->histOfs[1] = (icellX1*ncells.height + icellY0)*nbins;
data->histWeights[1] = cellX*(1.f - cellY);
//下面計算出的結果為9
data->histOfs[2] = (icellX0*ncells.height + icellY1)*nbins;
data->histWeights[2] = (1.f - cellX)*cellY;
//下面計算出的結果為27
data->histOfs[3] = (icellX1*ncells.height + icellY1)*nbins;
data->histWeights[3] = cellX*cellY;
}
else
//滿足這個else條件說明icellY0取-1或者1,也就是說block縱座標在(0, 3.5)
//和(11.5, 15)之間.
//此時的畫素點對相鄰的2個cell有權重貢獻
{
data = &pixData[rawBlockSize + (count2++)];
if( (unsigned)icellY0 < (unsigned)ncells.height )
{
//(unsigned)-1等於127>2,所以此處滿足if條件時icellY0==1;
//icellY1==1;
icellY1 = icellY0;
cellY = 1.f - cellY;
}
//不滿足if條件時,icellY0==-1;icellY1==0;
//當然了,這2種情況下icellX0==0;icellX1==1;
data->histOfs[0] = (icellX0*ncells.height + icellY1)*nbins;
data->histWeights[0] = (1.f - cellX)*cellY;
data->histOfs[1] = (icellX1*ncells.height + icellY1)*nbins;
data->histWeights[1] = cellX*cellY;
data->histOfs[2] = data->histOfs[3] = 0;
data->histWeights[2] = data->histWeights[3] = 0;
}
}
//當block中橫座標滿足在(0, 3.5)和(11.5, 15)範圍內時,即
//icellX0==-1或==1
else
{
if( (unsigned)icellX0 < (unsigned)ncells.width )
{
//icellX1=icllX0=1;
icellX1 = icellX0;
cellX = 1.f - cellX;
}
//當icllY0=0時,此時對2個cell有貢獻
if( (unsigned)icellY0 < (unsigned)ncells.height &&
(unsigned)icellY1 < (unsigned)ncells.height )
{
data = &pixData[rawBlockSize + (count2++)];
data->histOfs[0] = (icellX1*ncells.height + icellY0)*nbins;
data->histWeights[0] = cellX*(1.f - cellY);
data->histOfs[1] = (icellX1*ncells.height + icellY1)*nbins;
data->histWeights[1] = cellX*cellY;
data->histOfs[2] = data->histOfs[3] = 0;
data->histWeights[2] = data->histWeights[3] = 0;
}
else
//此時只對自身的cell有貢獻
{
data = &pixData[count1++];
if( (unsigned)icellY0 < (unsigned)ncells.height )
{
icellY1 = icellY0;
cellY = 1.f - cellY;
}
data->histOfs[0] = (icellX1*ncells.height + icellY1)*nbins;
data->histWeights[0] = cellX*cellY;
data->histOfs[1] = data->histOfs[2] = data->histOfs[3] = 0;
data->histWeights[1] = data->histWeights[2] = data->histWeights[3] = 0;
}
}
//為什麼每個block中i,j位置的gradOfs和qangleOfs都相同且是如下的計算公式呢?
//那是因為輸入的_img引數不是代表整幅圖片而是檢測視窗大小的圖片,所以每個
//檢測視窗中關於block的資訊可以看做是相同的
data->gradOfs = (grad.cols*i + j)*2;
data->qangleOfs = (qangle.cols*i + j)*2;
//每個block中i,j位置的權重都是固定的
data->gradWeight = weights(i,j);
}
//保證所有的點都被掃描了一遍
assert( count1 + count2 + count4 == rawBlockSize );
// defragment pixData
//將pixData中按照記憶體排滿,這樣節省了2/3的記憶體
for( j = 0; j < count2; j++ )
pixData[j + count1] = pixData[j + rawBlockSize];
for( j = 0; j < count4; j++ )
pixData[j + count1 + count2] = pixData[j + rawBlockSize*2];
//此時count2表示至多對2個cell有貢獻的所有畫素點的個數
count2 += count1;
//此時count4表示至多對4個cell有貢獻的所有畫素點的個數
count4 += count2;
//上面是初始化pixData,下面開始初始化blockData
// initialize blockData
for( j = 0; j < nblocks.width; j++ )
for( i = 0; i < nblocks.height; i++ )
{
BlockData& data = blockData[j*nblocks.height + i];
//histOfs表示該block對檢測視窗貢獻的hog描述變數起點在整個
//變數中的座標
data.histOfs = (j*nblocks.height + i)*blockHistogramSize;
//imgOffset表示該block的左上角在檢測視窗中的座標
data.imgOffset = Point(j*blockStride.width,i*blockStride.height);
}
//一個檢測視窗對應一個blockData記憶體,一個block對應一個pixData記憶體。
}
//pt為該block左上角在滑動視窗中的座標,buf為指向檢測視窗中blocData的指標
//函式返回一個block描述子的指標
const float* HOGCache::getBlock(Point pt, float* buf)
{
float* blockHist = buf;
assert(descriptor != 0);
Size blockSize = descriptor->blockSize;
pt += imgoffset;
CV_Assert( (unsigned)pt.x <= (unsigned)(grad.cols - blockSize.width) &&
(unsigned)pt.y <= (unsigned)(grad.rows - blockSize.height) );
if( useCache )
{
//cacheStride可以認為和blockStride是一樣的
//保證所獲取到HOGCache是我們所需要的,即在block移動過程中會出現
CV_Assert( pt.x % cacheStride.width == 0 &&
pt.y % cacheStride.height == 0 );
//cacheIdx表示的是block個數的座標
Point cacheIdx(pt.x/cacheStride.width,
(pt.y/cacheStride.height) % blockCache.rows);
//ymaxCached的長度為一個檢測視窗垂直方向上容納的block個數
if( pt.y != ymaxCached[cacheIdx.y] )
{
//取出blockCacheFlags的第cacheIdx.y行並且賦值為0
Mat_<uchar> cacheRow = blockCacheFlags.row(cacheIdx.y);
cacheRow = (uchar)0;
ymaxCached[cacheIdx.y] = pt.y;
}
//blockHist指向該點對應block所貢獻的hog描述子向量,初始值為空
blockHist = &blockCache[cacheIdx.y][cacheIdx.x*blockHistogramSize];
uchar& computedFlag = blockCacheFlags(cacheIdx.y, cacheIdx.x);
if( computedFlag != 0 )
return blockHist;
computedFlag = (uchar)1; // set it at once, before actual computing
}
int k, C1 = count1, C2 = count2, C4 = count4;
//
const float* gradPtr = (const float*)(grad.data + grad.step*pt.y) + pt.x*2;
const uchar* qanglePtr = qangle.data + qangle.step*pt.y + pt.x*2;
CV_Assert( blockHist != 0 );
#ifdef HAVE_IPP
ippsZero_32f(blockHist,blockHistogramSize);
#else
for( k = 0; k < blockHistogramSize; k++ )
blockHist[k] = 0.f;
#endif
const PixData* _pixData = &pixData[0];
//C1表示只對自己所在cell有貢獻的點的個數
for( k = 0; k < C1; k++ )
{
const PixData& pk = _pixData[k];
//a表示的是幅度指標
const float* a = gradPtr + pk.gradOfs;
float w = pk.gradWeight*pk.histWeights[0];
//h表示的是相位指標
const uchar* h = qanglePtr + pk.qangleOfs;
//幅度有2個通道是因為每個畫素點的幅值被分解到了其相鄰的兩個bin上了
//相位有2個通道是因為每個畫素點的相位的相鄰處都有的2個bin的序號
int h0 = h[0], h1 = h[1];
float* hist = blockHist + pk.histOfs[0];
float t0 = hist[h0] + a[0]*w;
float t1 = hist[h1] + a[1]*w;
//hist中放的為加權的梯度值
hist[h0] = t0; hist[h1] = t1;
}
for( ; k < C2; k++ )
{
const PixData& pk = _pixData[k];
const float* a = gradPtr + pk.gradOfs;
float w, t0, t1, a0 = a[0], a1 = a[1];
const uchar* h = qanglePtr + pk.qangleOfs;
int h0 = h[0], h1 = h[1];
//因為此時的畫素對2個cell有貢獻,這是其中一個cell的貢獻
float* hist = blockHist + pk.histOfs[0];
w = pk.gradWeight*pk.histWeights[0];
t0 = hist[h0] + a0*w;
t1 = hist[h1] + a1*w;
hist[h0] = t0; hist[h1] = t1;
//另一個cell的貢獻
hist = blockHist + pk.histOfs[1];
w = pk.gradWeight*pk.histWeights[1];
t0 = hist[h0] + a0*w;
t1 = hist[h1] + a1*w;
hist[h0] = t0; hist[h1] = t1;
}
//和上面類似
for( ; k < C4; k++ )
{
const PixData& pk = _pixData[k];
const float* a = gradPtr + pk.gradOfs;
float w, t0, t1, a0 = a[0], a1 = a[1];
const uchar* h = qanglePtr + pk.qangleOfs;
int h0 = h[0], h1 = h[1];
float* hist = blockHist + pk.histOfs[0];
w = pk.gradWeight*pk.histWeights[0];
t0 = hist[h0] + a0*w;
t1 = hist[h1] + a1*w;
hist[h0] = t0; hist[h1] = t1;
hist = blockHist + pk.histOfs[1];
w = pk.gradWeight*pk.histWeights[1];
t0 = hist[h0] + a0*w;
t1 = hist[h1] + a1*w;
hist[h0] = t0; hist[h1] = t1;
hist = blockHist + pk.histOfs[2];
w = pk.gradWeight*pk.histWeights[2];
t0 = hist[h0] + a0*w;
t1 = hist[h1] + a1*w;
hist[h0] = t0; hist[h1] = t1;
hist = blockHist + pk.histOfs[3];
w = pk.gradWeight*pk.histWeights[3];
t0 = hist[h0] + a0*w;
t1 = hist[h1] + a1*w;
hist[h0] = t0; hist[h1] = t1;
}
normalizeBlockHistogram(blockHist);
return blockHist;
}
void HOGCache::normalizeBlockHistogram(float* _hist) const
{
float* hist = &_hist[0];
#ifdef HAVE_IPP
size_t sz = blockHistogramSize;
#else
size_t i, sz = blockHistogramSize;
#endif
float sum = 0;
#ifdef HAVE_IPP
ippsDotProd_32f(hist,hist,sz,&sum);
#else
//第一次歸一化求的是平方和
for( i = 0; i < sz; i++ )
sum += hist[i]*hist[i];
#endif
//分母為平方和開根號+0.1
float scale = 1.f/(std::sqrt(sum)+sz*0.1f), thresh = (float)descriptor->L2HysThreshold;
#ifdef HAVE_IPP
ippsMulC_32f_I(scale,hist,sz);
ippsThreshold_32f_I( hist, sz, thresh, ippCmpGreater );
ippsDotProd_32f(hist,hist,sz,&sum);
#else
for( i = 0, sum = 0; i < sz; i++ )
{
//第2次歸一化是在第1次的基礎上繼續求平和和
hist[i] = std::min(hist[i]*scale, thresh);
sum += hist[i]*hist[i];
}
#endif
scale = 1.f/(std::sqrt(sum)+1e-3f);
#ifdef HAVE_IPP
ippsMulC_32f_I(scale,hist,sz);
#else
//最終歸一化結果
for( i = 0; i < sz; i++ )
hist[i] *= scale;
#endif
}
//返回測試圖片中水平方向和垂直方向共有多少個檢測視窗
Size HOGCache::windowsInImage(Size imageSize, Size winStride) const
{
return Size((imageSize.width - winSize.width)/winStride.width + 1,
(imageSize.height - winSize.height)/winStride.height + 1);
}
//給定圖片的大小,已經檢測視窗滑動的大小和測試圖片中的檢測視窗的索引,得到該索引處
//檢測視窗的尺寸,包括座標資訊
Rect HOGCache::getWindow(Size imageSize, Size winStride, int idx) const
{
int nwindowsX = (imageSize.width - winSize.width)/winStride.width + 1;
int y = idx / nwindowsX;//商
int x = idx - nwindowsX*y;//餘數
return Rect( x*winStride.width, y*winStride.height, winSize.width, winSize.height );
}
void HOGDescriptor::compute(const Mat& img, vector<float>& descriptors,
Size winStride, Size padding,
const vector<Point>& locations) const
{
//Size()表示長和寬都是0
if( winStride == Size() )
winStride = cellSize;
//gcd為求最大公約數,如果採用預設值的話,則2者相同
Size cacheStride(gcd(winStride.width, blockStride.width),
gcd(winStride.height, blockStride.height));
size_t nwindows = locations.size();
//alignSize(m, n)返回n的倍數大於等於m的最小值
padding.width = (int)alignSize(std::max(padding.width, 0), cacheStride.width);
padding.height = (int)alignSize(std::max(padding.height, 0), cacheStride.height);
Size paddedImgSize(img.cols + padding.width*2, img.rows + padding.height*2);
HOGCache cache(this, img, padding, padding, nwindows == 0, cacheStride);
if( !nwindows )
//Mat::area()表示為Mat的面積
nwindows = cache.windowsInImage(paddedImgSize, winStride).area();
const HOGCache::BlockData* blockData = &cache.blockData[0];
int nblocks = cache.nblocks.area();
int blockHistogramSize = cache.blockHistogramSize;
size_t dsize = getDescriptorSize();//一個hog的描述長度
//resize()為改變矩陣的行數,如果減少矩陣的行數則只保留減少後的
//那些行,如果是增加行數,則保留所有的行。
//這裡將描述子長度擴充套件到整幅圖片
descriptors.resize(dsize*nwindows);
for( size_t i = 0; i < nwindows; i++ )
{
//descriptor為第i個檢測視窗的描述子首位置。
float* descriptor = &descriptors[i*dsize];
Point pt0;
//非空
if( !locations.empty() )
{
pt0 = locations[i];
//非法的點
if( pt0.x < -padding.width || pt0.x > img.cols + padding.width - winSize.width ||
pt0.y < -padding.height || pt0.y > img.rows + padding.height - winSize.height )
continue;
}
//locations為空
else
{
//pt0為沒有擴充前影象對應的第i個檢測視窗
pt0 = cache.getWindow(paddedImgSize, winStride, (int)i).tl() - Point(padding);
CV_Assert(pt0.x % cacheStride.width == 0 && pt0.y % cacheStride.height == 0);
}
for( int j = 0; j < nblocks; j++ )
{
const HOGCache::BlockData& bj = blockData[j];
//pt為block的左上角相對檢測圖片的座標
Point pt = pt0 + bj.imgOffset;
//dst為該block在整個測試圖片的描述子的位置
float* dst = descriptor + bj.histOfs;
const float* src = cache.getBlock(pt, dst);
if( src != dst )
#ifdef HAVE_IPP
ippsCopy_32f(src,dst,blockHistogramSize);
#else
for( int k = 0; k < blockHistogramSize; k++ )
dst[k] = src[k];
#endif
}
}
}
void HOGDescriptor::detect(const Mat& img,
vector<Point>& hits, vector<double>& weights, double hitThreshold,
Size winStride, Size padding, const vector<Point>& locations) const
{
//hits裡面存的是符合檢測到目標的視窗的左上角頂點座標
hits.clear();
if( svmDetector.empty() )
return;
if( winStride == Size() )
winStride = cellSize;
Size cacheStride(gcd(winStride.width, blockStride.width),
gcd(winStride.height, blockStride.height));
size_t nwindows = locations.size();
padding.width = (int)alignSize(std::max(padding.width, 0), cacheStride.width);
padding.height = (int)alignSize(std::max(padding.height, 0), cacheStride.height);
Size paddedImgSize(img.cols + padding.width*2, img.rows + padding.height*2);
HOGCache cache(this, img, padding, padding, nwindows == 0, cacheStride);
if( !nwindows )
nwindows = cache.windowsInImage(paddedImgSize, winStride).area();
const HOGCache::BlockData* blockData = &cache.blockData[0];
int nblocks = cache.nblocks.area();
int blockHistogramSize = cache.blockHistogramSize;
size_t dsize = getDescriptorSize();
double rho = svmDetector.size() > dsize ? svmDetector[dsize] : 0;
vector<float> blockHist(blockHistogramSize);
for( size_t i = 0; i < nwindows; i++ )
{
Point pt0;
if( !locations.empty() )
{
pt0 = locations[i];
if( pt0.x < -padding.width || pt0.x > img.cols + padding.width - winSize.width ||
pt0.y < -padding.height || pt0.y > img.rows + padding.height - winSize.height )
continue;
}
else
{
pt0 = cache.getWindow(paddedImgSize, winStride, (int)i).tl() - Point(padding);
CV_Assert(pt0.x % cacheStride.width == 0 && pt0.y % cacheStride.height == 0);
}
double s = rho;
//svmVec指向svmDetector最前面那個元素
const float* svmVec = &svmDetector[0];
#ifdef HAVE_IPP
int j;
#else
int j, k;
#endif
for( j = 0; j < nblocks; j++, svmVec += blockHistogramSize )
{
const HOGCache::BlockData& bj = blockData[j];
Point pt = pt0 + bj.imgOffset;
//vec為測試圖片pt處的block貢獻的描述子指標
const float* vec = cache.getBlock(pt, &blockHist[0]);
#ifdef HAVE_IPP
Ipp32f partSum;
ippsDotProd_32f(vec,svmVec,blockHistogramSize,&partSum);
s += (double)partSum;
#else
for( k = 0; k <= blockHistogramSize - 4; k += 4 )
//const float* svmVec = &svmDetector[0];
s += vec[k]*svmVec[k] + vec[k+1]*svmVec[k+1] +
vec[k+2]*svmVec[k+2] + vec[k+3]*svmVec[k+3];
for( ; k < blockHistogramSize; k++ )
s += vec[k]*svmVec[k];
#endif
}
if( s >= hitThreshold )
{
hits.push_back(pt0);
weights.push_back(s);
}
}
}
//不用保留檢測到目標的可信度,即權重
void HOGDescriptor::detect(const Mat& img, vector<Point>& hits, double hitThreshold,
Size winStride, Size padding, const vector<Point>& locations) const
{
vector<double> weightsV;
detect(img, hits, weightsV, hitThreshold, winStride, padding, locations);
}
struct HOGInvoker
{
HOGInvoker( const HOGDescriptor* _hog, const Mat& _img,
double _hitThreshold, Size _winStride, Size _padding,
const double* _levelScale, ConcurrentRectVector* _vec,
ConcurrentDoubleVector* _weights=0, ConcurrentDoubleVector* _scales=0 )
{
hog = _hog;
img = _img;
hitThreshold = _hitThreshold;
winStride = _winStride;
padding = _padding;
levelScale = _levelScale;
vec = _vec;
weights = _weights;
scales = _scales;
}
void operator()( const BlockedRange& range ) const
{
int i, i1 = range.begin(), i2 = range.end();
double minScale = i1 > 0 ? levelScale[i1] : i2 > 1 ? levelScale[i1+1] : std::max(img.cols, img.rows);
//將原圖片進行縮放
Size maxSz(cvCeil(img.cols/minScale), cvCeil(img.rows/minScale));
Mat smallerImgBuf(maxSz, img.type());
vector<Point> locations;
vector<double> hitsWeights;
for( i = i1; i < i2; i++ )
{
double scale = levelScale[i];
Size sz(cvRound(img.cols/scale), cvRound(img.rows/scale));
//smallerImg只是構造一個指標,並沒有複製資料
Mat smallerImg(sz, img.type(), smallerImgBuf.data);
//沒有尺寸縮放
if( sz == img.size() )
smallerImg = Mat(sz, img.type(), img.data, img.step);
//有尺寸縮放
else
resize(img, smallerImg, sz);
//該函式實際上是將返回的值存在locations和histWeights中
//其中locations存的是目標區域的左上角座標
hog->detect(smallerImg, locations, hitsWeights, hitThreshold, winStride, padding);
Size scaledWinSize = Size(cvRound(hog->winSize.width*scale), cvRound(hog->winSize.height*scale));
for( size_t j = 0; j < locations.size(); j++ )
{
//儲存目標區域
vec->push_back(Rect(cvRound(locations[j].x*scale),
cvRound(locations[j].y*scale),
scaledWinSize.width, scaledWinSize.height));
//儲存縮放尺寸
if (scales) {
scales->push_back(scale);
}
}
//儲存svm計算後的結果值
if (weights && (!hitsWeights.empty()))
{
for (size_t j = 0; j < locations.size(); j++)
{
weights->push_back(hitsWeights[j]);
}
}
}
}
const HOGDescriptor* hog;
Mat img;
double hitThreshold;
Size winStride;
Size padding;
const double* levelScale;
//typedef tbb::concurrent_vector<Rect> ConcurrentRectVector;
ConcurrentRectVector* vec;
//typedef tbb::concurrent_vector<double> ConcurrentDoubleVector;
ConcurrentDoubleVector* weights;
ConcurrentDoubleVector* scales;
};
void HOGDescriptor::detectMultiScale(
const Mat& img, vector<Rect>& foundLocations, vector<double>& foundWeights,