尺度空间理论与图像金字塔（二）
SIFT簡(jiǎn)介
整理一下方便閱讀，作者寫(xiě)的東西摘自論文，在此感謝xiaowei等的貢獻(xiàn)
?
- DoG尺度空間構(gòu)造（Scale-space extrema detection）http://blog.csdn.net/xiaowei_cqu/article/details/8067881
- 關(guān)鍵點(diǎn)搜索與定位（Keypoint localization）http://blog.csdn.net/xiaowei_cqu/article/details/8087239
- 方向賦值（Orientation assignment）http://blog.csdn.net/xiaowei_cqu/article/details/8096072
- 關(guān)鍵點(diǎn)描述（Keypoint descriptor）http://blog.csdn.net/xiaowei_cqu/article/details/8113565
- OpenCV實(shí)現(xiàn)：特征檢測(cè)器FeatureDetectorhttp://blog.csdn.net/xiaowei_cqu/article/details/8652096
- SIFT中LoG和DoG的比較http://blog.csdn.net/xiaowei_cqu/article/details/27692123
DoG尺度空間構(gòu)造
?
自然界中的物體隨著觀測(cè)尺度不同有不同的表現(xiàn)形態(tài)。例如我們形容建筑物用“米”，觀測(cè)分子、原子等用“納米”。更形象的例子比如Google地圖，滑動(dòng)鼠標(biāo)輪可以改變觀測(cè)地圖的尺度，看到的地圖繪制也不同；還有電影中的拉伸鏡頭等等……
尺度空間中各尺度圖像的模糊程度逐漸變大，能夠模擬人在距離目標(biāo)由近到遠(yuǎn)時(shí)目標(biāo)在視網(wǎng)膜上的形成過(guò)程。
尺度越大圖像越模糊。?
為什么要討論尺度空間？
用機(jī)器視覺(jué)系統(tǒng)分析未知場(chǎng)景時(shí)，計(jì)算機(jī)并不預(yù)先知道圖像中物體的尺度。我們需要同時(shí)考慮圖像在多尺度下的描述，獲知感興趣物體的最佳尺度。另外如果不同的尺度下都有同樣的關(guān)鍵點(diǎn)，那么在不同的尺度的輸入圖像下就都可以檢測(cè)出來(lái)關(guān)鍵點(diǎn)匹配，也就是尺度不變性。
圖像的尺度空間表達(dá)就是圖像在所有尺度下的描述。
?
尺度空間表達(dá)與金字塔多分辨率表達(dá)
?
高斯模糊
高斯核是唯一可以產(chǎn)生多尺度空間的核（《Scale-space theory: A basic tool for analysing structures at different scales》）。一個(gè)圖像的尺度空間L(x,y,σ) ,定義為原始圖像I(x,y)與一個(gè)可變尺度的2維高斯函數(shù)G(x,y,σ)卷積運(yùn)算。?
二維空間高斯函數(shù)：
尺度空間：
尺度是自然客觀存在的，不是主觀創(chuàng)造的。高斯卷積只是表現(xiàn)尺度空間的一種形式。
二維空間高斯函數(shù)是等高線從中心成正太分布的同心圓：
?
?
分布不為零的點(diǎn)組成卷積陣與原始圖像做變換，即每個(gè)像素值是周?chē)噜徬袼刂档母咚蛊骄Ｒ粋€(gè)5*5的高斯模版如下所示：
?
?
高斯模版是圓對(duì)稱(chēng)的，且卷積的結(jié)果使原始像素值有最大的權(quán)重，距離中心越遠(yuǎn)的相鄰像素值權(quán)重也越小。
在實(shí)際應(yīng)用中，在計(jì)算高斯函數(shù)的離散近似時(shí)，在大概3σ距離之外的像素都可以看作不起作用，這些像素的計(jì)算也就可以忽略。所以，通常程序只計(jì)算(6σ+1)*(6σ+1)就可以保證相關(guān)像素影響。高斯模糊另一個(gè)很厲害的性質(zhì)就是線性可分：使用二維矩陣變換的高斯模糊可以通過(guò)在水平和豎直方向各進(jìn)行一維高斯矩陣變換相加得到。
O(N^2*m*n)次乘法就縮減成了O(N*m*n)+O(N*m*n)次乘法。（N為高斯核大小，m,n為二維圖像高和寬）
其實(shí)高斯這一部分只需要簡(jiǎn)單了解就可以了，在OpenCV也只需要一句代碼：
?
[cpp]?view plain?copy
- GaussianBlur(dbl,?dbl,?Size(),?sig_diff,?sig_diff);??
?
金字塔多分辨率
金字塔是早期圖像多尺度的表示形式。圖像金字塔化一般包括兩個(gè)步驟：使用低通濾波器平滑圖像；對(duì)平滑圖像進(jìn)行降采樣（通常是水平，豎直方向1/2），從而得到一系列尺寸縮小的圖像。
上圖中（a）是對(duì)原始信號(hào)進(jìn)行低通濾波，（b）是降采樣得到的信號(hào)。
而對(duì)于二維圖像，一個(gè)傳統(tǒng)的金字塔中，每一層圖像由上一層分辨率的長(zhǎng)、寬各一半，也就是四分之一的像素組成：
多尺度和多分辨率
尺度空間表達(dá)和金字塔多分辨率表達(dá)之間最大的不同是：
- 尺度空間表達(dá)是由不同高斯核平滑卷積得到，在所有尺度上有相同的分辨率；
- 而金字塔多分辨率表達(dá)每層分辨率減少固定比率。
所以，金字塔多分辨率生成較快，且占用存儲(chǔ)空間少；而多尺度表達(dá)隨著尺度參數(shù)的增加冗余信息也變多。
多尺度表達(dá)的優(yōu)點(diǎn)在于圖像的局部特征可以用簡(jiǎn)單的形式在不同尺度上描述；而金字塔表達(dá)沒(méi)有理論基礎(chǔ)，難以分析圖像局部特征。
?
DoG（Difference of Gaussian）
?
高斯拉普拉斯LoG金字塔
結(jié)合尺度空間表達(dá)和金字塔多分辨率表達(dá)，就是在使用尺度空間時(shí)使用金字塔表示，也就是計(jì)算機(jī)視覺(jué)中最有名的拉普拉斯金子塔（《The Laplacian pyramid as a compact image code》）。
高斯拉普拉斯LoG（Laplace of Guassian）算子就是對(duì)高斯函數(shù)進(jìn)行拉普拉斯變換：
核心思想還是高斯，這個(gè)不多敘述。
?
高斯差分DoG金字塔
DoG（Difference of Gaussian）其實(shí)是對(duì)高斯拉普拉斯LoG的近似，也就是對(duì)的近似。SIFT算法建議，在某一尺度上的特征檢測(cè)可以通過(guò)對(duì)兩個(gè)相鄰高斯尺度空間的圖像相減，得到DoG的響應(yīng)值圖像D(x,y,σ)。然后仿照LoG方法，通過(guò)對(duì)響應(yīng)值圖像D(x,y,σ)進(jìn)行局部最大值搜索，在空間位置和尺度空間定位局部特征點(diǎn)。其中：
k為相鄰兩個(gè)尺度空間倍數(shù)的常數(shù)。
上圖中（a）是DoG的三維圖，（b）是DoG與LoG的對(duì)比。
?
金字塔構(gòu)建
?
構(gòu)建高斯金字塔
為了得到DoG圖像，先要構(gòu)造高斯金字塔。我們回過(guò)頭來(lái)繼續(xù)說(shuō)高斯金字塔~
高斯金字塔在多分辨率金字塔簡(jiǎn)單降采樣基礎(chǔ)上加了高斯濾波，也就是對(duì)金字塔每層圖像用不同參數(shù)的σ做高斯模糊，使得每層金字塔有多張高斯模糊圖像。金字塔每層多張圖像合稱(chēng)為一組（Octave），每組有多張（也叫層Interval）圖像。另外，降采樣時(shí)，金字塔上邊一組圖像的第一張圖像（最底層的一張）是由前一組（金字塔下面一組）圖像的倒數(shù)第三張隔點(diǎn)采樣得到。
?
以下是OpenCV中構(gòu)建高斯金字塔的代碼，我加了相應(yīng)的注釋：
[cpp]?view plain?copy
- //?構(gòu)建nOctaves組（每組nOctaves+3層）高斯金字塔??
- void?SIFT::buildGaussianPyramid(?const?Mat&?base,?vector<Mat>&?pyr,?int?nOctaves?)?const??
- {??
- ????vector<double>?sig(nOctaveLayers?+?3);??
- ????pyr.resize(nOctaves*(nOctaveLayers?+?3));??
- ??
- ????//?precompute?Gaussian?sigmas?using?the?following?formula:??
- ????//??\sigma_{total}^2?=?\sigma_{i}^2?+?\sigma_{i-1}^2、??
- ????//?計(jì)算對(duì)圖像做不同尺度高斯模糊的尺度因子??
- ????sig[0]?=?sigma;??
- ????double?k?=?pow(?2.,?1.?/?nOctaveLayers?);??
- ????for(?int?i?=?1;?i?<?nOctaveLayers?+?3;?i++?)??
- ????{??
- ????????double?sig_prev?=?pow(k,?(double)(i-1))*sigma;??
- ????????double?sig_total?=?sig_prev*k;??
- ????????sig[i]?=?std::sqrt(sig_total*sig_total?-?sig_prev*sig_prev);??
- ????}??
- ??
- ????for(?int?o?=?0;?o?<?nOctaves;?o++?)??
- ????{??
- ????????//?DoG金子塔需要nOctaveLayers+2層圖像來(lái)檢測(cè)nOctaves層尺度??
- ????????//?所以高斯金字塔需要nOctaveLayers+3層圖像得到nOctaveLayers+2層DoG金字塔??
- ????????for(?int?i?=?0;?i?<?nOctaveLayers?+?3;?i++?)??
- ????????{??
- ????????????//?dst為第o組（Octave）金字塔??
- ????????????Mat&?dst?=?pyr[o*(nOctaveLayers?+?3)?+?i];??
- ????????????//?第0組第0層為原始圖像??
- ????????????if(?o?==?0??&&??i?==?0?)??
- ????????????????dst?=?base;??
- ??????????????
- ????????????//?base?of?new?octave?is?halved?image?from?end?of?previous?octave??
- ????????????//?每一組第0副圖像時(shí)上一組倒數(shù)第三幅圖像隔點(diǎn)采樣得到??
- ????????????else?if(?i?==?0?)??
- ????????????{??
- ????????????????const?Mat&?src?=?pyr[(o-1)*(nOctaveLayers?+?3)?+?nOctaveLayers];??
- ????????????????resize(src,?dst,?Size(src.cols/2,?src.rows/2),??
- ???????????????????????0,?0,?INTER_NEAREST);??
- ????????????}??
- ????????????//?每一組第i副圖像是由第i-1副圖像進(jìn)行sig[i]的高斯模糊得到??
- ????????????//?也就是本組圖像在sig[i]的尺度空間下的圖像??
- ????????????else??
- ????????????{??
- ????????????????const?Mat&?src?=?pyr[o*(nOctaveLayers?+?3)?+?i-1];??
- ????????????????GaussianBlur(src,?dst,?Size(),?sig[i],?sig[i]);??
- ????????????}??
- ????????}??
- ????}??
- }??
高斯金字塔的組數(shù)為：
?代碼10-17行是計(jì)算高斯模糊的系數(shù)σ，具體關(guān)系如下：
其中，σ為尺度空間坐標(biāo)，s為每組中層坐標(biāo)，σ0為初始尺度，S為每組層數(shù)（一般為3~5）。根據(jù)這個(gè)公式，我們可以得到金字塔組內(nèi)各層尺度以及組間各圖像尺度關(guān)系。
組內(nèi)相鄰圖像尺度關(guān)系：
相鄰組間尺度關(guān)系：
所以，相鄰兩組的同一層尺度為2倍的關(guān)系。
最終尺度序列總結(jié)為：
o為金字塔組數(shù)，n為每組金字塔層數(shù)。
?
構(gòu)建DoG金字塔
構(gòu)建高斯金字塔之后，就是用金字塔相鄰圖像相減構(gòu)造DoG金字塔。
?
?下面為構(gòu)造DoG的代碼：
[cpp]?view plain?copy
- //?構(gòu)建nOctaves組（每組nOctaves+2層）高斯差分金字塔??
- void?SIFT::buildDoGPyramid(?const?vector<Mat>&?gpyr,?vector<Mat>&?dogpyr?)?const??
- {??
- ????int?nOctaves?=?(int)gpyr.size()/(nOctaveLayers?+?3);??
- ????dogpyr.resize(?nOctaves*(nOctaveLayers?+?2)?);??
- ??
- ????for(?int?o?=?0;?o?<?nOctaves;?o++?)??
- ????{??
- ????????for(?int?i?=?0;?i?<?nOctaveLayers?+?2;?i++?)??
- ????????{??
- ????????????//?第o組第i副圖像為高斯金字塔中第o組第i+1和i組圖像相減得到??
- ????????????const?Mat&?src1?=?gpyr[o*(nOctaveLayers?+?3)?+?i];??
- ????????????const?Mat&?src2?=?gpyr[o*(nOctaveLayers?+?3)?+?i?+?1];??
- ????????????Mat&?dst?=?dogpyr[o*(nOctaveLayers?+?2)?+?i];??
- ????????????subtract(src2,?src1,?dst,?noArray(),?CV_16S);??
- ????????}??
- ????}??
- }??
-
尺度空間理論
-
由前一步《DoG尺度空間構(gòu)造》，我們得到了DoG高斯差分金字塔：
如上圖的金字塔，高斯尺度空間金字塔中每組有五層不同尺度圖像，相鄰兩層相減得到四層DoG結(jié)果。關(guān)鍵點(diǎn)搜索就在這四層DoG圖像上尋找局部極值點(diǎn)。
DoG局部極值點(diǎn)
尋找DoG極值點(diǎn)時(shí)，每一個(gè)像素點(diǎn)和它所有的相鄰點(diǎn)比較，當(dāng)其大于（或小于）它的圖像域和尺度域的所有相鄰點(diǎn)時(shí)，即為極值點(diǎn)。如下圖所示，比較的范圍是個(gè)3×3的立方體：中間的檢測(cè)點(diǎn)和它同尺度的8個(gè)相鄰點(diǎn)，以及和上下相鄰尺度對(duì)應(yīng)的9×2個(gè)點(diǎn)——共26個(gè)點(diǎn)比較，以確保在尺度空間和二維圖像空間都檢測(cè)到極值點(diǎn)。?
在一組中，搜索從每組的第二層開(kāi)始，以第二層為當(dāng)前層，第一層和第三層分別作為立方體的的上下層；搜索完成后再以第三層為當(dāng)前層做同樣的搜索。所以每層的點(diǎn)搜索兩次。通常我們將組Octaves索引以-1開(kāi)始，則在比較時(shí)犧牲了-1組的第0層和第N組的最高層
高斯金字塔，DoG圖像及極值計(jì)算的相互關(guān)系如上圖所示。
?
二 ??關(guān)鍵點(diǎn)精確定位
以上極值點(diǎn)的搜索是在離散空間進(jìn)行搜索的，由下圖可以看到，在離散空間找到的極值點(diǎn)不一定是真正意義上的極值點(diǎn)。可以通過(guò)對(duì)尺度空間DoG函數(shù)進(jìn)行曲線擬合尋找極值點(diǎn)來(lái)減小這種誤差。
利用DoG函數(shù)在尺度空間的Taylor展開(kāi)式：
則極值點(diǎn)為：
程序中還除去了極值小于0.04的點(diǎn)。如下所示：
[cpp]?view plain?copy
- //?Detects?features?at?extrema?in?DoG?scale?space.??Bad?features?are?discarded??
- //?based?on?contrast?and?ratio?of?principal?curvatures.??
- //?在DoG尺度空間尋特征點(diǎn)（極值點(diǎn)）??
- void?SIFT::findScaleSpaceExtrema(?const?vector<Mat>&?gauss_pyr,?const?vector<Mat>&?dog_pyr,??
- ??????????????????????????????????vector<KeyPoint>&?keypoints?)?const??
- {??
- ????int?nOctaves?=?(int)gauss_pyr.size()/(nOctaveLayers?+?3);??
- ??????
- ????//?The?contrast?threshold?used?to?filter?out?weak?features?in?semi-uniform??
- ????//?(low-contrast)?regions.?The?larger?the?threshold,?the?less?features?are?produced?by?the?detector.??
- ????//?過(guò)濾掉弱特征的閾值?contrastThreshold默認(rèn)為0.04??
- ????int?threshold?=?cvFloor(0.5?*?contrastThreshold?/?nOctaveLayers?*?255?*?SIFT_FIXPT_SCALE);??
- ????const?int?n?=?SIFT_ORI_HIST_BINS;?//36??
- ????float?hist[n];??
- ????KeyPoint?kpt;??
- ??
- ????keypoints.clear();??
- ??
- ????for(?int?o?=?0;?o?<?nOctaves;?o++?)??
- ????????for(?int?i?=?1;?i?<=?nOctaveLayers;?i++?)??
- ????????{??
- ????????????int?idx?=?o*(nOctaveLayers+2)+i;??
- ????????????const?Mat&?img?=?dog_pyr[idx];??
- ????????????const?Mat&?prev?=?dog_pyr[idx-1];??
- ????????????const?Mat&?next?=?dog_pyr[idx+1];??
- ????????????int?step?=?(int)img.step1();??
- ????????????int?rows?=?img.rows,?cols?=?img.cols;??
- ??
- ????????????for(?int?r?=?SIFT_IMG_BORDER;?r?<?rows-SIFT_IMG_BORDER;?r++)??
- ????????????{??
- ????????????????const?short*?currptr?=?img.ptr<short>(r);??
- ????????????????const?short*?prevptr?=?prev.ptr<short>(r);??
- ????????????????const?short*?nextptr?=?next.ptr<short>(r);??
- ??
- ????????????????for(?int?c?=?SIFT_IMG_BORDER;?c?<?cols-SIFT_IMG_BORDER;?c++)??
- ????????????????{??
- ????????????????????int?val?=?currptr[c];??
- ??
- ????????????????????//?find?local?extrema?with?pixel?accuracy??
- ????????????????????//?尋找局部極值點(diǎn)，DoG中每個(gè)點(diǎn)與其所在的立方體周?chē)?6個(gè)點(diǎn)比較??
- ????????????????????//?if?（val比所有都大?或者?val比所有都小）??
- ????????????????????if(?std::abs(val)?>?threshold?&&??
- ???????????????????????((val?>?0?&&?val?>=?currptr[c-1]?&&?val?>=?currptr[c+1]?&&??
- ?????????????????????????val?>=?currptr[c-step-1]?&&?val?>=?currptr[c-step]?&&???
- ?????????????????????????val?>=?currptr[c-step+1]?&&?val?>=?currptr[c+step-1]?&&???
- ?????????????????????????val?>=?currptr[c+step]?&&?val?>=?currptr[c+step+1]?&&??
- ?????????????????????????val?>=?nextptr[c]?&&?val?>=?nextptr[c-1]?&&???
- ?????????????????????????val?>=?nextptr[c+1]?&&?val?>=?nextptr[c-step-1]?&&???
- ?????????????????????????val?>=?nextptr[c-step]?&&?val?>=?nextptr[c-step+1]?&&???
- ?????????????????????????val?>=?nextptr[c+step-1]?&&?val?>=?nextptr[c+step]?&&???
- ?????????????????????????val?>=?nextptr[c+step+1]?&&?val?>=?prevptr[c]?&&???
- ?????????????????????????val?>=?prevptr[c-1]?&&?val?>=?prevptr[c+1]?&&??
- ?????????????????????????val?>=?prevptr[c-step-1]?&&?val?>=?prevptr[c-step]?&&???
- ?????????????????????????val?>=?prevptr[c-step+1]?&&?val?>=?prevptr[c+step-1]?&&???
- ?????????????????????????val?>=?prevptr[c+step]?&&?val?>=?prevptr[c+step+1])?||??
- ????????????????????????(val?<?0?&&?val?<=?currptr[c-1]?&&?val?<=?currptr[c+1]?&&??
- ?????????????????????????val?<=?currptr[c-step-1]?&&?val?<=?currptr[c-step]?&&???
- ?????????????????????????val?<=?currptr[c-step+1]?&&?val?<=?currptr[c+step-1]?&&???
- ?????????????????????????val?<=?currptr[c+step]?&&?val?<=?currptr[c+step+1]?&&??
- ?????????????????????????val?<=?nextptr[c]?&&?val?<=?nextptr[c-1]?&&???
- ?????????????????????????val?<=?nextptr[c+1]?&&?val?<=?nextptr[c-step-1]?&&???
- ?????????????????????????val?<=?nextptr[c-step]?&&?val?<=?nextptr[c-step+1]?&&???
- ?????????????????????????val?<=?nextptr[c+step-1]?&&?val?<=?nextptr[c+step]?&&???
- ?????????????????????????val?<=?nextptr[c+step+1]?&&?val?<=?prevptr[c]?&&???
- ?????????????????????????val?<=?prevptr[c-1]?&&?val?<=?prevptr[c+1]?&&??
- ?????????????????????????val?<=?prevptr[c-step-1]?&&?val?<=?prevptr[c-step]?&&???
- ?????????????????????????val?<=?prevptr[c-step+1]?&&?val?<=?prevptr[c+step-1]?&&???
- ?????????????????????????val?<=?prevptr[c+step]?&&?val?<=?prevptr[c+step+1])))??
- ????????????????????{??
- ????????????????????????int?r1?=?r,?c1?=?c,?layer?=?i;??
- ??????????????????????????
- ????????????????????????//?關(guān)鍵點(diǎn)精確定位??
- ????????????????????????if(?!adjustLocalExtrema(dog_pyr,?kpt,?o,?layer,?r1,?c1,??
- ????????????????????????????????????????????????nOctaveLayers,?(float)contrastThreshold,??
- ????????????????????????????????????????????????(float)edgeThreshold,?(float)sigma)?)??
- ????????????????????????????continue;??
- ??????????????????????????
- ????????????????????????float?scl_octv?=?kpt.size*0.5f/(1?<<?o);??
- ????????????????????????//?計(jì)算梯度直方圖??
- ????????????????????????float?omax?=?calcOrientationHist(??
- ????????????????????????????gauss_pyr[o*(nOctaveLayers+3)?+?layer],??
- ????????????????????????????Point(c1,?r1),??
- ????????????????????????????cvRound(SIFT_ORI_RADIUS?*?scl_octv),??
- ????????????????????????????SIFT_ORI_SIG_FCTR?*?scl_octv,??
- ????????????????????????????hist,?n);??
- ????????????????????????float?mag_thr?=?(float)(omax?*?SIFT_ORI_PEAK_RATIO);??
- ????????????????????????for(?int?j?=?0;?j?<?n;?j++?)??
- ????????????????????????{??
- ????????????????????????????int?l?=?j?>?0???j?-?1?:?n?-?1;??
- ????????????????????????????int?r2?=?j?<?n-1???j?+?1?:?0;??
- ??
- ????????????????????????????if(?hist[j]?>?hist[l]??&&??hist[j]?>?hist[r2]??&&??hist[j]?>=?mag_thr?)??
- ????????????????????????????{??
- ????????????????????????????????float?bin?=?j?+?0.5f?*?(hist[l]-hist[r2])?/???
- ????????????????????????????????(hist[l]?-?2*hist[j]?+?hist[r2]);??
- ????????????????????????????????bin?=?bin?<?0???n?+?bin?:?bin?>=?n???bin?-?n?:?bin;??
- ????????????????????????????????kpt.angle?=?(float)((360.f/n)?*?bin);??
- ????????????????????????????????keypoints.push_back(kpt);??
- ????????????????????????????}??
- ????????????????????????}??
- ????????????????????}??
- ????????????????}??
- ????????????}??
- ????????}??
- }??
刪除邊緣效應(yīng)
除了DoG響應(yīng)較低的點(diǎn)，還有一些響應(yīng)較強(qiáng)的點(diǎn)也不是穩(wěn)定的特征點(diǎn)。DoG對(duì)圖像中的邊緣有較強(qiáng)的響應(yīng)值，所以落在圖像邊緣的點(diǎn)也不是穩(wěn)定的特征點(diǎn)。
一個(gè)平坦的DoG響應(yīng)峰值在橫跨邊緣的地方有較大的主曲率，而在垂直邊緣的地方有較小的主曲率。主曲率可以通過(guò)2×2的Hessian矩陣H求出：
D值可以通過(guò)求臨近點(diǎn)差分得到。H的特征值與D的主曲率成正比，具體可參見(jiàn)Harris角點(diǎn)檢測(cè)算法。
為了避免求具體的值，我們可以通過(guò)H將特征值的比例表示出來(lái)。令為最大特征值，為最小特征值，那么：
?
Tr(H)表示矩陣H的跡，Det(H)表示H的行列式。
令表示最大特征值與最小特征值的比值，則有：
上式與兩個(gè)特征值的比例有關(guān)。隨著主曲率比值的增加，也會(huì)增加。我們只需要去掉比率大于一定值的特征點(diǎn)。Lowe論文中去掉r=10的點(diǎn)。
[cpp]?view plain?copy
- //?Interpolates?a?scale-space?extremum's?location?and?scale?to?subpixel??
- //?accuracy?to?form?an?image?feature.??Rejects?features?with?low?contrast.??
- //?Based?on?Section?4?of?Lowe's?paper.??
- //?特征點(diǎn)精確定位??
- static?bool?adjustLocalExtrema(?const?vector<Mat>&?dog_pyr,?KeyPoint&?kpt,?int?octv,??
- ????????????????????????????????int&?layer,?int&?r,?int&?c,?int?nOctaveLayers,??
- ????????????????????????????????float?contrastThreshold,?float?edgeThreshold,?float?sigma?)??
- {??
- ????const?float?img_scale?=?1.f/(255*SIFT_FIXPT_SCALE);??
- ????const?float?deriv_scale?=?img_scale*0.5f;??
- ????const?float?second_deriv_scale?=?img_scale;??
- ????const?float?cross_deriv_scale?=?img_scale*0.25f;??
- ??
- ????float?xi=0,?xr=0,?xc=0,?contr;??
- ????int?i?=?0;??
- ??
- ????//三維子像元插值??
- ????for(?;?i?<?SIFT_MAX_INTERP_STEPS;?i++?)??
- ????{??
- ????????int?idx?=?octv*(nOctaveLayers+2)?+?layer;??
- ????????const?Mat&?img?=?dog_pyr[idx];??
- ????????const?Mat&?prev?=?dog_pyr[idx-1];??
- ????????const?Mat&?next?=?dog_pyr[idx+1];??
- ??
- ????????Vec3f?dD((img.at<short>(r,?c+1)?-?img.at<short>(r,?c-1))*deriv_scale,??
- ?????????????????(img.at<short>(r+1,?c)?-?img.at<short>(r-1,?c))*deriv_scale,??
- ?????????????????(next.at<short>(r,?c)?-?prev.at<short>(r,?c))*deriv_scale);??
- ??
- ????????float?v2?=?(float)img.at<short>(r,?c)*2;??
- ????????float?dxx?=?(img.at<short>(r,?c+1)?+???
- ????????????????img.at<short>(r,?c-1)?-?v2)*second_deriv_scale;??
- ????????float?dyy?=?(img.at<short>(r+1,?c)?+???
- ????????????????img.at<short>(r-1,?c)?-?v2)*second_deriv_scale;??
- ????????float?dss?=?(next.at<short>(r,?c)?+???
- ????????????????prev.at<short>(r,?c)?-?v2)*second_deriv_scale;??
- ????????float?dxy?=?(img.at<short>(r+1,?c+1)?-???
- ????????????????img.at<short>(r+1,?c-1)?-?img.at<short>(r-1,?c+1)?+???
- ????????????????img.at<short>(r-1,?c-1))*cross_deriv_scale;??
- ????????float?dxs?=?(next.at<short>(r,?c+1)?-???
- ????????????????next.at<short>(r,?c-1)?-?prev.at<short>(r,?c+1)?+???
- ????????????????prev.at<short>(r,?c-1))*cross_deriv_scale;??
- ????????float?dys?=?(next.at<short>(r+1,?c)?-???
- ????????????????next.at<short>(r-1,?c)?-?prev.at<short>(r+1,?c)?+???
- ????????????????prev.at<short>(r-1,?c))*cross_deriv_scale;??
- ??
- ????????Matx33f?H(dxx,?dxy,?dxs,??
- ??????????????????dxy,?dyy,?dys,??
- ??????????????????dxs,?dys,?dss);??
- ??
- ????????Vec3f?X?=?H.solve(dD,?DECOMP_LU);??
- ??
- ????????xi?=?-X[2];??
- ????????xr?=?-X[1];??
- ????????xc?=?-X[0];??
- ??
- ????????if(?std::abs(?xi?)?<?0.5f??&&??std::abs(?xr?)?<?0.5f??&&??std::abs(?xc?)?<?0.5f?)??
- ????????????break;??
- ??
- ????????//將找到的極值點(diǎn)對(duì)應(yīng)成像素（整數(shù)）??
- ????????c?+=?cvRound(?xc?);??
- ????????r?+=?cvRound(?xr?);??
- ????????layer?+=?cvRound(?xi?);??
- ??
- ????????if(?layer?<?1?||?layer?>?nOctaveLayers?||??
- ???????????c?<?SIFT_IMG_BORDER?||?c?>=?img.cols?-?SIFT_IMG_BORDER??||??
- ???????????r?<?SIFT_IMG_BORDER?||?r?>=?img.rows?-?SIFT_IMG_BORDER?)??
- ????????????return?false;??
- ????}??
- ??
- ????/*?ensure?convergence?of?interpolation?*/??
- ????//?SIFT_MAX_INTERP_STEPS:插值最大步數(shù)，避免插值不收斂，程序中默認(rèn)為5??
- ????if(?i?>=?SIFT_MAX_INTERP_STEPS?)??
- ????????return?false;??
- ??
- ????{??
- ????????int?idx?=?octv*(nOctaveLayers+2)?+?layer;??
- ????????const?Mat&?img?=?dog_pyr[idx];??
- ????????const?Mat&?prev?=?dog_pyr[idx-1];??
- ????????const?Mat&?next?=?dog_pyr[idx+1];??
- ????????Matx31f?dD((img.at<short>(r,?c+1)?-?img.at<short>(r,?c-1))*deriv_scale,??
- ???????????????????(img.at<short>(r+1,?c)?-?img.at<short>(r-1,?c))*deriv_scale,??
- ???????????????????(next.at<short>(r,?c)?-?prev.at<short>(r,?c))*deriv_scale);??
- ????????float?t?=?dD.dot(Matx31f(xc,?xr,?xi));??
- ??
- ????????contr?=?img.at<short>(r,?c)*img_scale?+?t?*?0.5f;??
- ????????if(?std::abs(?contr?)?*?nOctaveLayers?<?contrastThreshold?)??
- ????????????return?false;??
- ??
- ????????/*?principal?curvatures?are?computed?using?the?trace?and?det?of?Hessian?*/??
- ???????//利用Hessian矩陣的跡和行列式計(jì)算主曲率的比值??
- ???????float?v2?=?img.at<short>(r,?c)*2.f;??
- ????????float?dxx?=?(img.at<short>(r,?c+1)?+???
- ????????????????img.at<short>(r,?c-1)?-?v2)*second_deriv_scale;??
- ????????float?dyy?=?(img.at<short>(r+1,?c)?+???
- ????????????????img.at<short>(r-1,?c)?-?v2)*second_deriv_scale;??
- ????????float?dxy?=?(img.at<short>(r+1,?c+1)?-???
- ????????????????img.at<short>(r+1,?c-1)?-?img.at<short>(r-1,?c+1)?+???
- ????????????????img.at<short>(r-1,?c-1))?*?cross_deriv_scale;??
- ????????float?tr?=?dxx?+?dyy;??
- ????????float?det?=?dxx?*?dyy?-?dxy?*?dxy;??
- ??
- ????????//這里edgeThreshold可以在調(diào)用SIFT()時(shí)輸入；??
- ????????//其實(shí)代碼中定義了?static?const?float?SIFT_CURV_THR?=?10.f?可以直接使用??
- ????????if(?det?<=?0?||?tr*tr*edgeThreshold?>=?(edgeThreshold?+?1)*(edgeThreshold?+?1)*det?)??
- ????????????return?false;??
- ????}??
- ??
- ????kpt.pt.x?=?(c?+?xc)?*?(1?<<?octv);??
- ????kpt.pt.y?=?(r?+?xr)?*?(1?<<?octv);??
- ????kpt.octave?=?octv?+?(layer?<<?8)?+?(cvRound((xi?+?0.5)*255)?<<?16);??
- ????kpt.size?=?sigma*powf(2.f,?(layer?+?xi)?/?nOctaveLayers)*(1?<<?octv)*2;??
- ??
- ????return?true;??
- }??
三 ?方向賦值
OpenCV】SIFT原理與源碼分析：方向賦值
?
由前一篇《關(guān)鍵點(diǎn)搜索與定位》，我們已經(jīng)找到了關(guān)鍵點(diǎn)。為了實(shí)現(xiàn)圖像旋轉(zhuǎn)不變性，需要根據(jù)檢測(cè)到的關(guān)鍵點(diǎn)局部圖像結(jié)構(gòu)為特征點(diǎn)方向賦值。也就是在findScaleSpaceExtrema()函數(shù)里看到的alcOrientationHist()語(yǔ)句：
[cpp]?view plain?copy
- //?計(jì)算梯度直方圖??
- float?omax?=?calcOrientationHist(gauss_pyr[o*(nOctaveLayers+3)?+?layer],??
- ????????????????????????????????????????????????Point(c1,?r1),??
- ????????????????????????????????????????????????cvRound(SIFT_ORI_RADIUS?*?scl_octv),??
- ????????????????????????????????????????????????SIFT_ORI_SIG_FCTR?*?scl_octv,??
- ????????????????????????????????????????????????hist,?n);??
我們使用圖像的梯度直方圖法求關(guān)鍵點(diǎn)局部結(jié)構(gòu)的穩(wěn)定方向。梯度方向和幅值
在前文中，精確定位關(guān)鍵點(diǎn)后也找到改特征點(diǎn)的尺度值σ，根據(jù)這一尺度值，得到最接近這一尺度值的高斯圖像：
使用有限差分，計(jì)算以關(guān)鍵點(diǎn)為中心，以3×1.5σ為半徑的區(qū)域內(nèi)圖像梯度的幅角和幅值，公式如下：
梯度直方圖
在完成關(guān)鍵點(diǎn)鄰域內(nèi)高斯圖像梯度計(jì)算后，使用直方圖統(tǒng)計(jì)鄰域內(nèi)像素對(duì)應(yīng)的梯度方向和幅值。
有關(guān)直方圖的基礎(chǔ)知識(shí)可以參考《數(shù)字圖像直方圖》，可以看做是離散點(diǎn)的概率表示形式。此處方向直方圖的核心是統(tǒng)計(jì)以關(guān)鍵點(diǎn)為原點(diǎn)，一定區(qū)域內(nèi)的圖像像素點(diǎn)對(duì)關(guān)鍵點(diǎn)方向生成所作的貢獻(xiàn)。
梯度方向直方圖的橫軸是梯度方向角，縱軸是剃度方向角對(duì)應(yīng)的梯度幅值累加值。梯度方向直方圖將0°~360°的范圍分為36個(gè)柱，每10°為一個(gè)柱。下圖是從高斯圖像上求取梯度，再由梯度得到梯度方向直方圖的例圖。
在計(jì)算直方圖時(shí)，每個(gè)加入直方圖的采樣點(diǎn)都使用圓形高斯函數(shù)函數(shù)進(jìn)行了加權(quán)處理，也就是進(jìn)行高斯平滑。這主要是因?yàn)镾IFT算法只考慮了尺度和旋轉(zhuǎn)不變形，沒(méi)有考慮仿射不變性。通過(guò)高斯平滑，可以使關(guān)鍵點(diǎn)附近的梯度幅值有較大權(quán)重，從而部分彌補(bǔ)沒(méi)考慮仿射不變形產(chǎn)生的特征點(diǎn)不穩(wěn)定。
通常離散的梯度直方圖要進(jìn)行插值擬合處理，以求取更精確的方向角度值。（這和《關(guān)鍵點(diǎn)搜索與定位》中插值的思路是一樣的）。
關(guān)鍵點(diǎn)方向
直方圖峰值代表該關(guān)鍵點(diǎn)處鄰域內(nèi)圖像梯度的主方向，也就是該關(guān)鍵點(diǎn)的主方向。在梯度方向直方圖中，當(dāng)存在另一個(gè)相當(dāng)于主峰值 ? ?80%能量的峰值時(shí)，則將這個(gè)方向認(rèn)為是該關(guān)鍵點(diǎn)的輔方向。所以一個(gè)關(guān)鍵點(diǎn)可能檢測(cè)得到多個(gè)方向，這可以增強(qiáng)匹配的魯棒性。Lowe的論文指出大概有15%關(guān)鍵點(diǎn)具有多方向，但這些點(diǎn)對(duì)匹配的穩(wěn)定性至為關(guān)鍵。
獲得圖像關(guān)鍵點(diǎn)主方向后，每個(gè)關(guān)鍵點(diǎn)有三個(gè)信息(x,y,σ,θ)：位置、尺度、方向。由此我們可以確定一個(gè)SIFT特征區(qū)域。通常使用一個(gè)帶箭頭的圓或直接使用箭頭表示SIFT區(qū)域的三個(gè)值：中心表示特征點(diǎn)位置，半徑表示關(guān)鍵點(diǎn)尺度（r=2.5σ）,箭頭表示主方向。具有多個(gè)方向的關(guān)鍵點(diǎn)可以復(fù)制成多份，然后將方向值分別賦給復(fù)制后的關(guān)鍵點(diǎn)。如下圖：
源碼
[cpp]?view plain?copy
- //?Computes?a?gradient?orientation?histogram?at?a?specified?pixel??
- //?計(jì)算特定點(diǎn)的梯度方向直方圖??
- static?float?calcOrientationHist(?const?Mat&?img,?Point?pt,?int?radius,??
- ??????????????????????????????????float?sigma,?float*?hist,?int?n?)??
- {??
- ????//len：2r+1也就是以r為半徑的圓（正方形）像素個(gè)數(shù)??
- ????int?i,?j,?k,?len?=?(radius*2+1)*(radius*2+1);??
- ??
- ????float?expf_scale?=?-1.f/(2.f?*?sigma?*?sigma);??
- ????AutoBuffer<float>?buf(len*4?+?n+4);??
- ????float?*X?=?buf,?*Y?=?X?+?len,?*Mag?=?X,?*Ori?=?Y?+?len,?*W?=?Ori?+?len;??
- ????float*?temphist?=?W?+?len?+?2;??
- ??
- ????for(?i?=?0;?i?<?n;?i++?)??
- ????????temphist[i]?=?0.f;??
- ??
- ????//?圖像梯度直方圖統(tǒng)計(jì)的像素范圍??
- ????for(?i?=?-radius,?k?=?0;?i?<=?radius;?i++?)??
- ????{??
- ????????int?y?=?pt.y?+?i;??
- ????????if(?y?<=?0?||?y?>=?img.rows?-?1?)??
- ????????????continue;??
- ????????for(?j?=?-radius;?j?<=?radius;?j++?)??
- ????????{??
- ????????????int?x?=?pt.x?+?j;??
- ????????????if(?x?<=?0?||?x?>=?img.cols?-?1?)??
- ????????????????continue;??
- ??
- ????????????float?dx?=?(float)(img.at<short>(y,?x+1)?-?img.at<short>(y,?x-1));??
- ????????????float?dy?=?(float)(img.at<short>(y-1,?x)?-?img.at<short>(y+1,?x));??
- ??
- ????????????X[k]?=?dx;?Y[k]?=?dy;?W[k]?=?(i*i?+?j*j)*expf_scale;??
- ????????????k++;??
- ????????}??
- ????}??
- ??
- ????len?=?k;??
- ??
- ????//?compute?gradient?values,?orientations?and?the?weights?over?the?pixel?neighborhood??
- ????exp(W,?W,?len);???
- ????fastAtan2(Y,?X,?Ori,?len,?true);???
- ????magnitude(X,?Y,?Mag,?len);???
- ??????
- ????//?計(jì)算直方圖的每個(gè)bin??
- ????for(?k?=?0;?k?<?len;?k++?)??
- ????{??
- ????????int?bin?=?cvRound((n/360.f)*Ori[k]);??
- ????????if(?bin?>=?n?)??
- ????????????bin?-=?n;??
- ????????if(?bin?<?0?)??
- ????????????bin?+=?n;??
- ????????temphist[bin]?+=?W[k]*Mag[k];??
- ????}??
- ??
- ????//?smooth?the?histogram??
- ????//?高斯平滑??
- ????temphist[-1]?=?temphist[n-1];??
- ????temphist[-2]?=?temphist[n-2];??
- ????temphist[n]?=?temphist[0];??
- ????temphist[n+1]?=?temphist[1];??
- ????for(?i?=?0;?i?<?n;?i++?)??
- ????{??
- ????????hist[i]?=?(temphist[i-2]?+?temphist[i+2])*(1.f/16.f)?+??
- ????????????(temphist[i-1]?+?temphist[i+1])*(4.f/16.f)?+??
- ????????????temphist[i]*(6.f/16.f);??
- ????}??
- ??????
- ????//?得到主方向??
- ????float?maxval?=?hist[0];??
- ????for(?i?=?1;?i?<?n;?i++?)??
- ????????maxval?=?std::max(maxval,?hist[i]);??
- ??
- ????return?maxval;??
- }??
- ?
四 【OpenCV】SIFT原理與源碼分析：關(guān)鍵點(diǎn)描述
《SIFT原理與源碼分析》系列文章索引：http://blog.csdn.net/xiaowei_cqu/article/details/8069548
?
由前一篇《方向賦值》，為找到的關(guān)鍵點(diǎn)即SIFT特征點(diǎn)賦了值，包含位置、尺度和方向的信息。接下來(lái)的步驟是關(guān)鍵點(diǎn)描述，即用用一組向量將這個(gè)關(guān)鍵點(diǎn)描述出來(lái)，這個(gè)描述子不但包括關(guān)鍵點(diǎn)，也包括關(guān)鍵點(diǎn)周?chē)鷮?duì)其有貢獻(xiàn)的像素點(diǎn)。用來(lái)作為目標(biāo)匹配的依據(jù)（所以描述子應(yīng)該有較高的獨(dú)特性，以保證匹配率），也可使關(guān)鍵點(diǎn)具有更多的不變特性，如光照變化、3D視點(diǎn)變化等。
SIFT描述子h(x,y,θ)是對(duì)關(guān)鍵點(diǎn)附近鄰域內(nèi)高斯圖像梯度統(tǒng)計(jì)的結(jié)果，是一個(gè)三維矩陣，但通常用一個(gè)矢量來(lái)表示。矢量通過(guò)對(duì)三維矩陣按一定規(guī)律排列得到。
描述子采樣區(qū)域
特征描述子與關(guān)鍵點(diǎn)所在尺度有關(guān)，因此對(duì)梯度的求取應(yīng)在特征點(diǎn)對(duì)應(yīng)的高斯圖像上進(jìn)行。將關(guān)鍵點(diǎn)附近劃分成d×d個(gè)子區(qū)域，每個(gè)子區(qū)域尺寸為mσ個(gè)像元（d=4，m=3，σ為尺特征點(diǎn)的尺度值）。考慮到實(shí)際計(jì)算時(shí)需要雙線性插值，故計(jì)算的圖像區(qū)域?yàn)閙σ(d+1)，再考慮旋轉(zhuǎn)，則實(shí)際計(jì)算的圖像區(qū)域?yàn)?#xff0c;如下圖所示：
源碼
[cpp]?view plain?copy
- ???Point?pt(cvRound(ptf.x),?cvRound(ptf.y));??
- //計(jì)算余弦，正弦，CV_PI/180:將角度值轉(zhuǎn)化為幅度值??
- ???float?cos_t?=?cosf(ori*(float)(CV_PI/180));??
- ???float?sin_t?=?sinf(ori*(float)(CV_PI/180));??
- ???float?bins_per_rad?=?n?/?360.f;??
- ???float?exp_scale?=?-1.f/(d?*?d?*?0.5f);?//d:SIFT_DESCR_WIDTH?4??????
- ???float?hist_width?=?SIFT_DESCR_SCL_FCTR?*?scl;??//?SIFT_DESCR_SCL_FCTR:?3???
- ???????????????????????????????????????????????//?scl:?size*0.5f??
- //?計(jì)算圖像區(qū)域半徑mσ(d+1)/2*sqrt(2)??
- //?1.4142135623730951f?為根號(hào)2??
- ???int?radius?=?cvRound(hist_width?*?1.4142135623730951f?*?(d?+?1)?*?0.5f);??
- ???cos_t?/=?hist_width;??
- ???sin_t?/=?hist_width;??
?
區(qū)域坐標(biāo)軸旋轉(zhuǎn)
為了保證特征矢量具有旋轉(zhuǎn)不變性，要以特征點(diǎn)為中心，在附近鄰域內(nèi)旋轉(zhuǎn)θ角，即旋轉(zhuǎn)為特征點(diǎn)的方向。
旋轉(zhuǎn)后區(qū)域內(nèi)采樣點(diǎn)新的坐標(biāo)為：
源碼
[cpp]?view plain?copy
- //計(jì)算采樣區(qū)域點(diǎn)坐標(biāo)旋轉(zhuǎn)??
- ????for(?i?=?-radius,?k?=?0;?i?<=?radius;?i++?)??
- ????????for(?j?=?-radius;?j?<=?radius;?j++?)??
- ????????{??
- ????????????/*?
- ?????????????Calculate?sample's?histogram?array?coords?rotated?relative?to?ori.?
- ?????????????Subtract?0.5?so?samples?that?fall?e.g.?in?the?center?of?row?1?(i.e.?
- ?????????????r_rot?=?1.5)?have?full?weight?placed?in?row?1?after?interpolation.?
- ?????????????*/??
- ????????????float?c_rot?=?j?*?cos_t?-?i?*?sin_t;??
- ????????????float?r_rot?=?j?*?sin_t?+?i?*?cos_t;??
- ????????????float?rbin?=?r_rot?+?d/2?-?0.5f;??
- ????????????float?cbin?=?c_rot?+?d/2?-?0.5f;??
- ????????????int?r?=?pt.y?+?i,?c?=?pt.x?+?j;??
- ??
- ????????????if(?rbin?>?-1?&&?rbin?<?d?&&?cbin?>?-1?&&?cbin?<?d?&&??
- ???????????????r?>?0?&&?r?<?rows?-?1?&&?c?>?0?&&?c?<?cols?-?1?)??
- ????????????{??
- ????????????????float?dx?=?(float)(img.at<short>(r,?c+1)?-?img.at<short>(r,?c-1));??
- ????????????????float?dy?=?(float)(img.at<short>(r-1,?c)?-?img.at<short>(r+1,?c));??
- ????????????????X[k]?=?dx;?Y[k]?=?dy;?RBin[k]?=?rbin;?CBin[k]?=?cbin;??
- ????????????????W[k]?=?(c_rot?*?c_rot?+?r_rot?*?r_rot)*exp_scale;??
- ????????????????k++;??
- ????????????}??
- ????????}??
計(jì)算采樣區(qū)域梯度直方圖
將旋轉(zhuǎn)后區(qū)域劃分為d×d個(gè)子區(qū)域（每個(gè)區(qū)域間隔為mσ像元），在子區(qū)域內(nèi)計(jì)算8個(gè)方向的梯度直方圖，繪制每個(gè)方向梯度方向的累加值，形成一個(gè)種子點(diǎn)。
與求主方向不同的是，此時(shí)，每個(gè)子區(qū)域梯度方向直方圖將0°~360°劃分為8個(gè)方向區(qū)間，每個(gè)區(qū)間為45°。即每個(gè)種子點(diǎn)有8個(gè)方向區(qū)間的梯度強(qiáng)度信息。由于存在d×d，即4×4個(gè)子區(qū)域，所以最終共有4×4×8=128個(gè)數(shù)據(jù)，形成128維SIFT特征矢量。
對(duì)特征矢量需要加權(quán)處理，加權(quán)采用mσd/2的標(biāo)準(zhǔn)高斯函數(shù)。為了除去光照變化影響，還有一步歸一化處理。
源碼
[cpp]?view plain?copy
- //計(jì)算梯度直方圖??
- ????for(?k?=?0;?k?<?len;?k++?)??
- ????{??
- ????????float?rbin?=?RBin[k],?cbin?=?CBin[k];??
- ????????float?obin?=?(Ori[k]?-?ori)*bins_per_rad;??
- ????????float?mag?=?Mag[k]*W[k];??
- ??
- ????????int?r0?=?cvFloor(?rbin?);??
- ????????int?c0?=?cvFloor(?cbin?);??
- ????????int?o0?=?cvFloor(?obin?);??
- ????????rbin?-=?r0;??
- ????????cbin?-=?c0;??
- ????????obin?-=?o0;??
- ??
- ????????//n為SIFT_DESCR_HIST_BINS：8，即將360°分為8個(gè)區(qū)間??
- ????????if(?o0?<?0?)??
- ????????????o0?+=?n;??
- ????????if(?o0?>=?n?)??
- ????????????o0?-=?n;??
- ??????????
- ??
- ????????//?histogram?update?using?tri-linear?interpolation??
- ????????//?雙線性插值??
- ????????float?v_r1?=?mag*rbin,?v_r0?=?mag?-?v_r1;??
- ????????float?v_rc11?=?v_r1*cbin,?v_rc10?=?v_r1?-?v_rc11;??
- ????????float?v_rc01?=?v_r0*cbin,?v_rc00?=?v_r0?-?v_rc01;??
- ????????float?v_rco111?=?v_rc11*obin,?v_rco110?=?v_rc11?-?v_rco111;??
- ????????float?v_rco101?=?v_rc10*obin,?v_rco100?=?v_rc10?-?v_rco101;??
- ????????float?v_rco011?=?v_rc01*obin,?v_rco010?=?v_rc01?-?v_rco011;??
- ????????float?v_rco001?=?v_rc00*obin,?v_rco000?=?v_rc00?-?v_rco001;??
- ??
- ????????int?idx?=?((r0+1)*(d+2)?+?c0+1)*(n+2)?+?o0;??
- ????????hist[idx]?+=?v_rco000;??
- ????????hist[idx+1]?+=?v_rco001;??
- ????????hist[idx+(n+2)]?+=?v_rco010;??
- ????????hist[idx+(n+3)]?+=?v_rco011;??
- ????????hist[idx+(d+2)*(n+2)]?+=?v_rco100;??
- ????????hist[idx+(d+2)*(n+2)+1]?+=?v_rco101;??
- ????????hist[idx+(d+3)*(n+2)]?+=?v_rco110;??
- ????????hist[idx+(d+3)*(n+2)+1]?+=?v_rco111;??
- ????}??
關(guān)鍵點(diǎn)描述源碼
[cpp]?view plain?copy
- //?SIFT關(guān)鍵點(diǎn)特征描述??
- //?SIFT描述子是關(guān)鍵點(diǎn)領(lǐng)域高斯圖像提取統(tǒng)計(jì)結(jié)果的一種表示??
- static?void?calcSIFTDescriptor(?const?Mat&?img,?Point2f?ptf,?float?ori,?float?scl,??
- ???????????????????????????????int?d,?int?n,?float*?dst?)??
- ?????????????????????????????
- {??
- ????Point?pt(cvRound(ptf.x),?cvRound(ptf.y));??
- ????//計(jì)算余弦，正弦，CV_PI/180:將角度值轉(zhuǎn)化為幅度值??
- ????float?cos_t?=?cosf(ori*(float)(CV_PI/180));??
- ????float?sin_t?=?sinf(ori*(float)(CV_PI/180));??
- ????float?bins_per_rad?=?n?/?360.f;??
- ????float?exp_scale?=?-1.f/(d?*?d?*?0.5f);?//d:SIFT_DESCR_WIDTH?4?????
- ????float?hist_width?=?SIFT_DESCR_SCL_FCTR?*?scl;??//?SIFT_DESCR_SCL_FCTR:?3???
- ???????????????????????????????????????????????????//?scl:?size*0.5f??
- ????//?計(jì)算圖像區(qū)域半徑mσ(d+1)/2*sqrt(2)??
- ????//?1.4142135623730951f?為根號(hào)2??
- ????int?radius?=?cvRound(hist_width?*?1.4142135623730951f?*?(d?+?1)?*?0.5f);??
- ????cos_t?/=?hist_width;??
- ????sin_t?/=?hist_width;??
- ??
- ????int?i,?j,?k,?len?=?(radius*2+1)*(radius*2+1),?histlen?=?(d+2)*(d+2)*(n+2);??
- ????int?rows?=?img.rows,?cols?=?img.cols;??
- ??
- ????AutoBuffer<float>?buf(len*6?+?histlen);??
- ????float?*X?=?buf,?*Y?=?X?+?len,?*Mag?=?Y,?*Ori?=?Mag?+?len,?*W?=?Ori?+?len;??
- ????float?*RBin?=?W?+?len,?*CBin?=?RBin?+?len,?*hist?=?CBin?+?len;??
- ??
- ????//初始化直方圖??
- ????for(?i?=?0;?i?<?d+2;?i++?)??
- ????{??
- ????????for(?j?=?0;?j?<?d+2;?j++?)??
- ????????????for(?k?=?0;?k?<?n+2;?k++?)??
- ????????????????hist[(i*(d+2)?+?j)*(n+2)?+?k]?=?0.;??
- ????}??
- ??
- ????//計(jì)算采樣區(qū)域點(diǎn)坐標(biāo)旋轉(zhuǎn)??
- ????for(?i?=?-radius,?k?=?0;?i?<=?radius;?i++?)??
- ????????for(?j?=?-radius;?j?<=?radius;?j++?)??
- ????????{??
- ????????????/*?
- ?????????????Calculate?sample's?histogram?array?coords?rotated?relative?to?ori.?
- ?????????????Subtract?0.5?so?samples?that?fall?e.g.?in?the?center?of?row?1?(i.e.?
- ?????????????r_rot?=?1.5)?have?full?weight?placed?in?row?1?after?interpolation.?
- ?????????????*/??
- ????????????float?c_rot?=?j?*?cos_t?-?i?*?sin_t;??
- ????????????float?r_rot?=?j?*?sin_t?+?i?*?cos_t;??
- ????????????float?rbin?=?r_rot?+?d/2?-?0.5f;??
- ????????????float?cbin?=?c_rot?+?d/2?-?0.5f;??
- ????????????int?r?=?pt.y?+?i,?c?=?pt.x?+?j;??
- ??
- ????????????if(?rbin?>?-1?&&?rbin?<?d?&&?cbin?>?-1?&&?cbin?<?d?&&??
- ???????????????r?>?0?&&?r?<?rows?-?1?&&?c?>?0?&&?c?<?cols?-?1?)??
- ????????????{??
- ????????????????float?dx?=?(float)(img.at<short>(r,?c+1)?-?img.at<short>(r,?c-1));??
- ????????????????float?dy?=?(float)(img.at<short>(r-1,?c)?-?img.at<short>(r+1,?c));??
- ????????????????X[k]?=?dx;?Y[k]?=?dy;?RBin[k]?=?rbin;?CBin[k]?=?cbin;??
- ????????????????W[k]?=?(c_rot?*?c_rot?+?r_rot?*?r_rot)*exp_scale;??
- ????????????????k++;??
- ????????????}??
- ????????}??
- ??
- ????len?=?k;??
- ????fastAtan2(Y,?X,?Ori,?len,?true);??
- ????magnitude(X,?Y,?Mag,?len);??
- ????exp(W,?W,?len);??
- ??
- ??????
- ????//計(jì)算梯度直方圖??
- ????for(?k?=?0;?k?<?len;?k++?)??
- ????{??
- ????????float?rbin?=?RBin[k],?cbin?=?CBin[k];??
- ????????float?obin?=?(Ori[k]?-?ori)*bins_per_rad;??
- ????????float?mag?=?Mag[k]*W[k];??
- ??
- ????????int?r0?=?cvFloor(?rbin?);??
- ????????int?c0?=?cvFloor(?cbin?);??
- ????????int?o0?=?cvFloor(?obin?);??
- ????????rbin?-=?r0;??
- ????????cbin?-=?c0;??
- ????????obin?-=?o0;??
- ??
- ????????//n為SIFT_DESCR_HIST_BINS：8，即將360°分為8個(gè)區(qū)間??
- ????????if(?o0?<?0?)??
- ????????????o0?+=?n;??
- ????????if(?o0?>=?n?)??
- ????????????o0?-=?n;??
- ??????????
- ??
- ????????//?histogram?update?using?tri-linear?interpolation??
- ????????//?雙線性插值??
- ????????float?v_r1?=?mag*rbin,?v_r0?=?mag?-?v_r1;??
- ????????float?v_rc11?=?v_r1*cbin,?v_rc10?=?v_r1?-?v_rc11;??
- ????????float?v_rc01?=?v_r0*cbin,?v_rc00?=?v_r0?-?v_rc01;??
- ????????float?v_rco111?=?v_rc11*obin,?v_rco110?=?v_rc11?-?v_rco111;??
- ????????float?v_rco101?=?v_rc10*obin,?v_rco100?=?v_rc10?-?v_rco101;??
- ????????float?v_rco011?=?v_rc01*obin,?v_rco010?=?v_rc01?-?v_rco011;??
- ????????float?v_rco001?=?v_rc00*obin,?v_rco000?=?v_rc00?-?v_rco001;??
- ??
- ????????int?idx?=?((r0+1)*(d+2)?+?c0+1)*(n+2)?+?o0;??
- ????????hist[idx]?+=?v_rco000;??
- ????????hist[idx+1]?+=?v_rco001;??
- ????????hist[idx+(n+2)]?+=?v_rco010;??
- ????????hist[idx+(n+3)]?+=?v_rco011;??
- ????????hist[idx+(d+2)*(n+2)]?+=?v_rco100;??
- ????????hist[idx+(d+2)*(n+2)+1]?+=?v_rco101;??
- ????????hist[idx+(d+3)*(n+2)]?+=?v_rco110;??
- ????????hist[idx+(d+3)*(n+2)+1]?+=?v_rco111;??
- ????}??
- ??
- ????//?finalize?histogram,?since?the?orientation?histograms?are?circular??
- ????//?最后確定直方圖，目標(biāo)方向直方圖是圓的??
- ????for(?i?=?0;?i?<?d;?i++?)??
- ????????for(?j?=?0;?j?<?d;?j++?)??
- ????????{??
- ????????????int?idx?=?((i+1)*(d+2)?+?(j+1))*(n+2);??
- ????????????hist[idx]?+=?hist[idx+n];??
- ????????????hist[idx+1]?+=?hist[idx+n+1];??
- ????????????for(?k?=?0;?k?<?n;?k++?)??
- ????????????????dst[(i*d?+?j)*n?+?k]?=?hist[idx+k];??
- ????????}??
- ????//?copy?histogram?to?the?descriptor,??
- ????//?apply?hysteresis?thresholding??
- ????//?and?scale?the?result,?so?that?it?can?be?easily?converted??
- ????//?to?byte?array??
- ????float?nrm2?=?0;??
- ????len?=?d*d*n;??
- ????for(?k?=?0;?k?<?len;?k++?)??
- ????????nrm2?+=?dst[k]*dst[k];??
- ????float?thr?=?std::sqrt(nrm2)*SIFT_DESCR_MAG_THR;??
- ????for(?i?=?0,?nrm2?=?0;?i?<?k;?i++?)??
- ????{??
- ????????float?val?=?std::min(dst[i],?thr);??
- ????????dst[i]?=?val;??
- ????????nrm2?+=?val*val;??
- ????}??
- ????nrm2?=?SIFT_INT_DESCR_FCTR/std::max(std::sqrt(nrm2),?FLT_EPSILON);??
- ????for(?k?=?0;?k?<?len;?k++?)??
- ????{??
- ????????dst[k]?=?saturate_cast<uchar>(dst[k]*nrm2);??
- ????}??
- } ?
五 【OpenCV】特征檢測(cè)器 FeatureDetector
OpenCV提供FeatureDetector實(shí)現(xiàn)特征檢測(cè)及匹配
?
[cpp]?view plain?copy
- class?CV_EXPORTS?FeatureDetector??
- {??
- public:??
- ????virtual?~FeatureDetector();??
- ????void?detect(?const?Mat&?image,?vector<KeyPoint>&?keypoints,??
- ????????const?Mat&?mask=Mat()?)?const;??
- ????void?detect(?const?vector<Mat>&?images,??
- ????????vector<vector<KeyPoint>?>&?keypoints,??
- ????????const?vector<Mat>&?masks=vector<Mat>()?)?const;??
- ????virtual?void?read(const?FileNode&);??
- ????virtual?void?write(FileStorage&)?const;??
- ????static?Ptr<FeatureDetector>?create(?const?string&?detectorType?);??
- protected:??
- ????...??
- };??
FeatureDetetor是虛類(lèi)，通過(guò)定義FeatureDetector的對(duì)象可以使用多種特征檢測(cè)方法。通過(guò)create()函數(shù)調(diào)用：?
?
[cpp]?view plain?copy
- Ptr<FeatureDetector>?FeatureDetector::create(const?string&?detectorType);??
?
OpenCV 2.4.3提供了10種特征檢測(cè)方法：
?
- "FAST" – FastFeatureDetector
- "STAR" – StarFeatureDetector
- "SIFT" – SIFT (nonfree module)
- "SURF" – SURF (nonfree module)
- "ORB" – ORB
- "MSER" – MSER
- "GFTT" – GoodFeaturesToTrackDetector
- "HARRIS" – GoodFeaturesToTrackDetector with Harris detector enabled
- "Dense" – DenseFeatureDetector
- "SimpleBlob" – SimpleBlobDetector
圖片中的特征大體可分為三種：點(diǎn)特征、線特征、塊特征。
FAST算法是Rosten提出的一種快速提取的點(diǎn)特征[1]，Harris與GFTT也是點(diǎn)特征，更具體來(lái)說(shuō)是角點(diǎn)特征（參考這里）。
SimpleBlob是簡(jiǎn)單塊特征，可以通過(guò)設(shè)置SimpleBlobDetector的參數(shù)決定提取圖像塊的主要性質(zhì)，提供5種：
顏色?By color、面積?By area、圓形度?By circularity、最大inertia?（不知道怎么翻譯）與最小inertia的比例?By ratio of the minimum inertia to maximum inertia、以及凸性?By convexity.
最常用的當(dāng)屬SIFT，尺度不變特征匹配算法（參考這里）；以及后來(lái)發(fā)展起來(lái)的SURF，都可以看做較為復(fù)雜的塊特征。這兩個(gè)算法在OpenCV nonfree的模塊里面，需要在附件引用項(xiàng)中添加opencv_nonfree243.lib，同時(shí)在代碼中加入：
[cpp]?view plain?copy
- initModule_nonfree();??
至于其他幾種算法，我就不太了解了?^_^
?
一個(gè)簡(jiǎn)單的使用演示：
[cpp]?view plain?copy
- int?main()??
- {??
- ??
- ????initModule_nonfree();//if?use?SIFT?or?SURF??
- ????Ptr<FeatureDetector>?detector?=?FeatureDetector::create(?"SIFT"?);??
- ????Ptr<DescriptorExtractor>?descriptor_extractor?=?DescriptorExtractor::create(?"SIFT"?);??
- ????Ptr<DescriptorMatcher>?descriptor_matcher?=?DescriptorMatcher::create(?"BruteForce"?);??
- ????if(?detector.empty()?||?descriptor_extractor.empty()?)??
- ????????throw?runtime_error("fail?to?create?detector!");??
- ??
- ????Mat?img1?=?imread("images\\box_in_scene.png");??
- ????Mat?img2?=?imread("images\\box.png");??
- ??
- ????//detect?keypoints;??
- ????vector<KeyPoint>?keypoints1,keypoints2;??
- ????detector->detect(?img1,?keypoints1?);??
- ????detector->detect(?img2,?keypoints2?);??
- ????cout?<<"img1:"<<?keypoints1.size()?<<?"?points??img2:"?<<keypoints2.size()???
- ????????<<?"?points"?<<?endl?<<?">"?<<?endl;??
- ??
- ????//compute?descriptors?for?keypoints;??
- ????cout?<<?"<?Computing?descriptors?for?keypoints?from?images..."?<<?endl;??
- ????Mat?descriptors1,descriptors2;??
- ????descriptor_extractor->compute(?img1,?keypoints1,?descriptors1?);??
- ????descriptor_extractor->compute(?img2,?keypoints2,?descriptors2?);??
- ??
- ????cout<<endl<<"Descriptors?Size:?"<<descriptors2.size()<<"?>"<<endl;??
- ????cout<<endl<<"Descriptor's?Column:?"<<descriptors2.cols<<endl??
- ????????<<"Descriptor's?Row:?"<<descriptors2.rows<<endl;??
- ????cout?<<?">"?<<?endl;??
- ??
- ????//Draw?And?Match?img1,img2?keypoints??
- ????Mat?img_keypoints1,img_keypoints2;??
- ????drawKeypoints(img1,keypoints1,img_keypoints1,Scalar::all(-1),0);??
- ????drawKeypoints(img2,keypoints2,img_keypoints2,Scalar::all(-1),0);??
- ????imshow("Box_in_scene?keyPoints",img_keypoints1);??
- ????imshow("Box?keyPoints",img_keypoints2);??
- ??
- ????descriptor_extractor->compute(?img1,?keypoints1,?descriptors1?);????
- ????vector<DMatch>?matches;??
- ????descriptor_matcher->match(?descriptors1,?descriptors2,?matches?);??
- ??
- ????Mat?img_matches;??
- ????drawMatches(img1,keypoints1,img2,keypoints2,matches,img_matches,Scalar::all(-1),CV_RGB(255,255,255),Mat(),4);??
- ??
- ????imshow("Mathc",img_matches);??
- ????waitKey(10000);??
- ????return?0;??
- }??
特征檢測(cè)結(jié)果如圖：Box_in_scene
?
Box
?
特征點(diǎn)匹配結(jié)果：
Match
?
另一點(diǎn)需要一提的是SimpleBlob的實(shí)現(xiàn)是有Bug的。不能直接通過(guò)?Ptr<FeatureDetector> detector = FeatureDetector::create("SimpleBlob"); ?語(yǔ)句來(lái)調(diào)用，而應(yīng)該直接創(chuàng)建?SimpleBlobDetector的對(duì)象：
[cpp]?view plain?copy
- ???????Mat?image?=?imread("images\\features.jpg");??
- Mat?descriptors;??
- vector<KeyPoint>?keypoints;??
- SimpleBlobDetector::Params?params;??
- //params.minThreshold?=?10;??
- //params.maxThreshold?=?100;??
- //params.thresholdStep?=?10;??
- //params.minArea?=?10;???
- //params.minConvexity?=?0.3;??
- //params.minInertiaRatio?=?0.01;??
- //params.maxArea?=?8000;??
- //params.maxConvexity?=?10;??
- //params.filterByColor?=?false;??
- //params.filterByCircularity?=?false;??
- SimpleBlobDetector?blobDetector(?params?);??
- blobDetector.create("SimpleBlob");??
- blobDetector.detect(?image,?keypoints?);??
- drawKeypoints(image,?keypoints,?image,?Scalar(255,0,0));??
六 計(jì)算機(jī)視覺(jué)】SIFT中LoG和DoG比較
實(shí)際計(jì)算時(shí)，三種方法計(jì)算的金字塔組數(shù)noctaves，尺度空間坐標(biāo)σ，以及每組金字塔內(nèi)的層數(shù)S是一樣的。同時(shí)，假設(shè)圖像為640*480的標(biāo)準(zhǔn)圖像。
金字塔層數(shù)：
其中o_min = 0，對(duì)于分辨率為640*480的圖像N=5。
每組金字塔內(nèi)圖像數(shù)：
S=3，即在做極值檢測(cè)時(shí)使用金子塔內(nèi)中間3張圖像。
對(duì)于LoG每組金字塔內(nèi)有S+2張圖像（S=-1，0,1,2,3），需要做S+1次高斯模糊操作（后一張圖像由前一張做高斯模糊得到）；而DoG每組金字塔有S+3張高斯圖像，得到S+2張DoG圖像。
尺度空間系數(shù)：其中，S表示每組金字塔內(nèi)圖像層數(shù)，n為當(dāng)前高斯層數(shù)，取0-4。DoG需要5個(gè)尺度系數(shù)得到6張GSS圖像，而LoG只需要前4個(gè)尺度系數(shù)得到5張圖像。
LoG
高斯核使用正太分布（高斯函數(shù)）計(jì)算模糊模版，N維空間正太分布方程為：
?
于是，二維高斯模板上的距離中心點(diǎn)為(x,y)的元素對(duì)應(yīng)高斯計(jì)算公式為：?
規(guī)范化的高斯拉普拉斯圖像為?
最終構(gòu)造LoG金字塔有5層，每層有S+2=5張圖像，每層金字塔內(nèi)每張圖像尺度是前一張的k倍，即構(gòu)成的連續(xù)尺度序列：?
其中o為當(dāng)前金字塔層數(shù)，n為在當(dāng)前金字塔層中圖像張數(shù)。
由于卷積計(jì)算性質(zhì)：在計(jì)算時(shí)，通過(guò)對(duì)前一張圖像做尺度系數(shù)為的卷積操作，可以減少卷積計(jì)算次數(shù)。故在金字塔每層內(nèi)的S+2張圖像，需要S+1次卷積操作，每次LoG核的尺度系數(shù)為：
?
圖1.?LoG金字塔示意圖。 左側(cè)為圖像尺度空間系數(shù)，紅色矩形框?yàn)閷?shí)際參與比較的空間系數(shù)，右側(cè)為在上一張圖像做LoG卷積操作的尺度系數(shù)
?
DoG
由于LoG與Gauss核具有如下關(guān)系：
?
即，
?
因此，LoG算子可以用高斯差分算子DoG（Difference of Guassians）表示。?
于是通過(guò)高斯金字塔每層內(nèi)相鄰兩張圖像相減可以得到DoG金字塔。對(duì)于最后需要S張圖像尋找極值點(diǎn)，需要S+2張DoG圖像，S+3張高斯圖像。具體關(guān)系如圖2.所示。?
圖 2. 由高斯金字塔得到DoG金字塔及其對(duì)應(yīng)的尺度空間系數(shù)示意圖。
?
LoG & DoG
一個(gè)直觀的比較結(jié)果，使用分別計(jì)算LoG和DoG，得到：
可以看到，LoG比DoG明顯需要更多的加法運(yùn)算和乘法運(yùn)算。雖然DoG需要在每層金字塔多做一次高斯操作（即為了得到S+2張DoG圖需要S+3張高斯模糊圖），但通過(guò)減法取代LoG核的計(jì)算過(guò)程，顯著減少了運(yùn)算次數(shù)，大大節(jié)省了運(yùn)算時(shí)間。
總結(jié)
以上是生活随笔為你收集整理的尺度空间理论与图像金字塔（二）的全部?jī)?nèi)容，希望文章能夠幫你解決所遇到的問(wèn)題。
- 上一篇： 肝脏移植手术多少钱？
- 下一篇： 求欧若拉的歌词