0 引言
特徵提取就是從影像中提取顯著並且具有可區分性和可匹配性的點結構。常見的點結構一般為影像內容中的角點、交叉點、閉合區域中心點等具有一定物理結構的點，而提取點結構的一般思想為構建能夠區分其他影像結構的響應函式或者從特徵線或輪廓中進行稀疏取樣。Harris角點檢測器便是運用二階矩或自相關矩陣來加速區域性極值搜尋並保證方向的不變性。基於畫素比較的特徵提取方法也稱為二值特徵，通常具有極高的提取效率並具有一定的方向不變性以及所提取的特徵點具有較高的重複率，對後續的匹配具有重要意義，然而這類方法受尺度和仿射變換的影響較大。針對上述問題，帶有尺度資訊的斑點特徵成為特徵提取的另一種形式，其最早是由Lindeberg 等人提出的高斯拉普拉斯（Laplace of Gaussian，LoG）函式響應來實現，並從中提出了尺度空間理論，其利用高斯響應函式的圓對稱性和對區域性團結構的極值響應特性以及對噪聲抑制能力，透過不同高斯標準差實現在尺度空間上的極值搜尋，從而提取對尺度、方向和噪聲魯棒的特徵點並得到相應的尺度資訊。為了避免大量的計算，D.Lowe 等人介紹了一種高斯差分（Difference-of-Gaussian，DoG）法來近似LoG的計算，並提出了著名的SIFT特徵描述子。
SIFT是由Lowe[2004]開發的一個演算法，用於提取影像中的不變特徵。稱它為變換的原因是，它會將影像資料變換為相對於區域性影像特徵的尺度不變座標。SIFT是本章到目前為止討論的最複雜的特徵檢測和描述方法。本節中用到了大量由實驗確定的引數。因此，與前面介紹的許多方法不同，SIFT具有很強的探索性，因為我們目前所具備的知識還不能讓我們將各種方法組合為一個能夠求解單重方法所不能求解的問題的“系統”。因此，我們不得不透過實驗確定控制更復雜系統效能的各種引數的相互作用。SIFT特徵（稱為關鍵點）對影像尺度和旋轉是不變的，並且對仿射失真、三維視點變化、噪聲和光照變化具有很強的魯棒性。SIFT的輸入是一幅影像，輸出是一個n維特徵向量，向量的元素是不變的特徵描述子。
1 SIFT演算法流程
Lowe將SIFT演算法分為如下四步：
（1）尺度空間極值檢測：搜尋所有尺度上的影像位置。透過高斯差分函式來識別潛在的對於尺度和旋轉不變的關鍵點。
（2）關鍵點定位：在每個候選的位置上，透過一個擬合精細的模型來確定位置和尺度。關鍵點的選擇依據於它們的穩定程度。
（3）關鍵點方向確定：基於影像區域性的梯度方向，分配給每個關鍵點位置一個或多個方向。所有後面的對影像資料的操作都相對於關鍵點的方向、尺度和位置進項變換，從而保證了對於這些變換的不變性。
（4）關鍵點描述：在每個關鍵點周圍的鄰域內，在選定的尺度上測量影像區域性的梯度。這些梯度作為關鍵點的描述符，它允許比較大的區域性形狀的變形或光照變化。
1.1 尺度空間極值檢測
1.1.1 尺度空間
SIFT演算法的第一階段是找到對尺度變化不變的影像位置。將影像利用不同的尺度引數進行平滑後得到一堆影像，目的是模擬影像的尺度減小出現的細節損失。在SIFT中，用於實現平滑的是高斯核，因此尺度引數是標準差。標準差的大小決定影像的平滑程度，大尺度對應影像的概貌特徵，小尺度對應影像的細節特徵。大的值對應粗糙尺度(低解析度)，小的值對應精細尺度(高解析度)。因此灰度影像f(x,y)的尺度空間 L(x,y,σ)是f與一個可變尺度高斯核G(x,y,σ)的卷積：
式中，尺度由引數σ控制，G的形式如下：
輸入影像f(x,y)依次與標準差為σ，kσ，k2σ，k3σ，...，的高斯核卷積，生成一堆由一個常量因子k分隔的高斯濾波（平滑）影像。
Q1: 怎麼理解平滑建立的尺度空間？
對於一副影像，近距離觀察和遠距離觀察看到的影像效果是不同的，前者比較清晰，透過前者能看到影像的一些細節資訊，後者比較模糊，透過後者能看到影像的一些輪廓的資訊，這就是影像的尺度，影像的尺度是自然存在的，並不是人為創造的。尺度空間中各尺度影像的模糊程度逐漸變大，能夠模擬人在距離目標由近到遠時目標在視網膜上的形成過程。在平滑後影像上的一點，包含的是原圖中一塊畫素的資訊，屬於語義層面的多尺度。
1.1.2 高斯金字塔
影像金字塔廣泛應用於各種視覺應用中。影像金字塔是影像的集合，它是由原始影像產生，連續降取樣，直到達到一些期望的停止點。常出現的兩種金字塔是：高斯金字塔和拉普拉斯金字塔。高斯金字塔用於降取樣，當我們要從金字塔中較低的影像重構上取樣影像時，需要拉普拉斯金字塔。SIFT中的高斯金字塔與最原始的高斯金字塔有些區別，因為它在構造尺度空間時，將這些不同尺度影像分為了O個Octave、每個Octave又分為了S層。
- 對於O的選擇：O = log2(min(M,N))-3，M、N指原影像的長和寬，求最小值後開log再減3
- 對於S的選擇：S = n + 3，其中，n指我們希望提取多少個圖片中的特徵。當S=2時，每個Octave有5層。為什麼要+3呢？因為後面要差分，5層差分得到4層，然後每層要和上下兩層比較，也就是4層中有兩層沒法比較，只留下兩層，即n是2。
- 關於比例因子k：為了滿足尺度變化的連續性，某一組的每一層應該相差一個k倍。對於Octave2中的第一層，我們取Octave1中的倒數第三層，因為倒數第三層的σ為knσ，也就是為了湊2σ，達到一個隔點取點的降取樣效果。所以有knσ=2σ，得k=1/n，當n=2時，k等於√2。即Octave1：【σ，...，knσ，kn+1σ，kn+2σ，kn+3σ】,Octave2：【knσ，kn+1σ，kn+2σ，kn+3σ,kn+4σ】
Q2:按照尺度空間的理論，模糊就已經改變尺度了，為什麼還要下采樣呢？
void SIFT_Impl::buildGaussianPyramid(const Mat& base, std::vector<Mat>& pyr, int nOctaves) const
{
CV_TRACE_FUNCTION();
std::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
sig[0] = sigma;
double k = std::pow(2., 1. / nOctaveLayers);
for (int i = 1; i < nOctaveLayers + 3; i++)
{
double sig_prev = std::pow(k, (double)(i - 1)) * sigma;
double sig_total = sig_prev * k;
//為了降低運算，將大尺度的高斯運算元拆分為2個小尺度高斯運算元的卷積
sig[i] = std::sqrt(sig_total * sig_total - sig_prev * sig_prev);
}
for (int o = 0; o < nOctaves; o++)
{
for (int i = 0; i < nOctaveLayers + 3; i++)
{
Mat& dst = pyr[o * (nOctaveLayers + 3) + i];
if (o == 0 && i == 0)
dst = base;
// base of new octave is halved image from end of previous octave
else if (i == 0)//每個Octave的第一層，都是上一個Octave的倒數第三層的下采樣，這樣可以得到連續的尺度
{
const Mat& src = pyr[(o - 1) * (nOctaveLayers + 3) + nOctaveLayers];
resize(src, dst, Size(src.cols / 2, src.rows / 2),
0, 0, INTER_NEAREST);
}
else
{ //每組中下一層由上一層高斯模糊得到，直接對第一層高斯模糊計算量太大
const Mat& src = pyr[o * (nOctaveLayers + 3) + i - 1];
GaussianBlur(src, dst, Size(), sig[i], sig[i]);
}
}
}
}
1.1.3 高斯拉普拉斯(LoG)金字塔
在介紹邊緣檢測的部落格中，我們瞭解到二階導數的過零點可以用於確定粗邊緣的中心位置。Marr-Hildreth邊緣檢測子利用高斯拉普拉斯（LoG）函式來檢測邊緣，檢測運算元可以為不同尺度，大運算元檢測模糊邊緣，小運算元檢測清晰的細節。LoG也是一種比較常用的斑點檢測方法。
上面一行原始訊號類似於一條線的灰度剖面，下面是在相同尺度σ下的LoG結果。可以看出總的LoG曲線其實是兩條邊界上LoG函式的結果的疊加，當兩條邊界足夠接近時，可以將這個原始訊號看作一個blob，兩個LoG函式的結果也合二為一得到一個極值點，這個極值點對應blob的中心。所以邊緣檢測對應的是LoG的過零點，而斑點檢測對應的是LoG的極值點。
檢測blob時，需要選擇合適的尺度，然後LOG可以透過檢測極值點來檢測blob。上圖第一行為原訊號對不同尺度LoG的響應，可以發現隨著尺度的不斷增大，LoG曲線由雙波谷逐漸融合成單波谷，但是響應的幅值越來越弱。這是因為，隨著尺度的增大，LoG運算元的最大幅度逐漸減小，導致響應也隨著尺度的增大而減小。因此應該使用σ2對LoG進行歸一化，即σ2▽2G。上圖第二行是歸一化後的拉普拉斯響應。歸一化後選擇產生最強響應的尺度，在該尺度上對應的極值就是blob的中心位置。理論表明，對於一個圓形blob，當二維LoG運算元的零點值曲線和目標圓形邊緣重合時取得最強響應。
即令LoG=0，
得到σ = r/√2。使用LoG尋找斑點時不僅在影像上尋找極值點，還要求在尺度空間上也是極值點。也就是檢測在尺度空間和影像空間都是極值的點，就是blob區域的中心點。
1.1.4 高斯差分(DOG)金字塔
二元連續高斯函式G(x,y,σ)對σ求導，可得
而LoG函式為
兩者僅相差一個σ，即
根據極限又可以得到
因此，當k≈1時，有
兩個不同方差的高斯函式相減記作DoG，和尺度歸一化的LoG只相差一個倍數k-1，這個倍數是一個常數，不影響極值點檢測。DoG計算速度快，LoG更精確。因此，我們用高斯差DoG近似高斯拉普拉斯LoG，具體做法就是將相鄰兩層的高斯尺度影像相減：
為了在每組中檢測n個尺度的極值點，DOG金字塔每組需n+2層影像，因為DOG金字塔的每一層和其上下兩層比較得到極值點，第一層和最後一層是不能檢測極值點的，而DOG金字塔由高斯金字塔相鄰兩層相減得到，則高斯金字塔每組需n+3層影像。
OpenCV構建DoG金字塔原始碼
平行計算，共nOctaves個Octave，每個Octave有nOctaveLayers+3層，共nOctaves * (nOctaveLayers + 2)個差分計算。
class buildDoGPyramidComputer : public ParallelLoopBody
{
public:
buildDoGPyramidComputer(
int _nOctaveLayers,
const std::vector<Mat>& _gpyr,
std::vector<Mat>& _dogpyr)
: nOctaveLayers(_nOctaveLayers),
gpyr(_gpyr),
dogpyr(_dogpyr) { }
void operator()(const cv::Range& range) const CV_OVERRIDE
{
CV_TRACE_FUNCTION();
const int begin = range.start;
const int end = range.end;
for (int a = begin; a < end; a++)
{
const int o = a / (nOctaveLayers + 2);//確定該影像屬於哪個Octave
const int i = a % (nOctaveLayers + 2);//確定該影像屬於Octave的第幾層
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(), DataType<sift_wt>::type);//第i+1層減去第i層
}
}
private:
int nOctaveLayers;
const std::vector<Mat>& gpyr;
std::vector<Mat>& dogpyr;
};
void SIFT_Impl::buildDoGPyramid(const std::vector<Mat>& gpyr, std::vector<Mat>& dogpyr) const
{
CV_TRACE_FUNCTION();
int nOctaves = (int)gpyr.size() / (nOctaveLayers + 3);
dogpyr.resize(nOctaves * (nOctaveLayers + 2));//儲存所有的DoG影像
// 平行計算，共nOctaves個Octave，每個Octave有nOctaveLayers+3層，共nOctaves * (nOctaveLayers + 2)個差分計算
parallel_for_(Range(0, nOctaves * (nOctaveLayers + 2)), buildDoGPyramidComputer(nOctaveLayers, gpyr, dogpyr));
}
1.2 找關鍵點
1.2.1 查詢初始關鍵點
由高斯金字塔得到每組n+3層影像，所有相鄰兩幅高斯濾波影像得到n+2個差函式D(x,y,σ)。將這些差函式視為影像，影像的細節水平隨著我們上調尺度空間而下降。下圖顯示了SIFT查詢影像D(x,y,σ)中的極值的過程。
在D(x,y,σ)影像中的每個位置（顯示為黑色x），將該位置的畫素值與其在當前影像中的8個相鄰畫素值及其在上方和下方影像中的9個相鄰畫素值，即26個相鄰畫素，進行比較。如果該位置的值大於其所有相鄰畫素的值，或小於所有相鄰畫素的值，那麼該位置被選為極值（最大值或最小值點）。在一個Octave的第一個（最後一個）尺度中檢測不到任何極值，因為它沒有相同大小的下方（上方）尺度影像。
OpenCV查詢初始關鍵點原始碼
在DoG尺度空間中找初始極值點，是透過比較每一個畫素與其26個鄰域畫素相比較。如果該位置的值大於其所有相鄰畫素的值，或小於所有相鄰畫素的值，那麼該位置被選為極值。OpenCV中對應的是findScaleSpaceExtremaT類，其中包括了極值計算，調整極值，梯度計算等。下面為極值計算部分。
void process(const cv::Range& range)
{
CV_TRACE_FUNCTION();
const int begin = range.start;
const int end = range.end;
static const int n = SIFT_ORI_HIST_BINS;
float CV_DECL_ALIGNED(CV_SIMD_WIDTH) hist[n];
const Mat& img = dog_pyr[idx];
const Mat& prev = dog_pyr[idx - 1];
const Mat& next = dog_pyr[idx + 1];
for (int r = begin; r < end; r++)
{
const sift_wt* currptr = img.ptr<sift_wt>(r);
const sift_wt* prevptr = prev.ptr<sift_wt>(r);
const sift_wt* nextptr = next.ptr<sift_wt>(r);
int c = SIFT_IMG_BORDER;
// vector loop reminder, better predictibility and less branch density
for (; c < cols - SIFT_IMG_BORDER; c++)
{
sift_wt val = currptr[c];
if (std::abs(val) <= threshold)
continue;
sift_wt _00, _01, _02;
sift_wt _10, _12;
sift_wt _20, _21, _22;
_00 = currptr[c - step - 1]; _01 = currptr[c - step]; _02 = currptr[c - step + 1];
_10 = currptr[c - 1]; _12 = currptr[c + 1];
_20 = currptr[c + step - 1]; _21 = currptr[c + step]; _22 = currptr[c + step + 1];
bool calculate = false;
if (val > 0)//當val大於0，找極大值
{
sift_wt vmax = std::max(std::max(std::max(_00, _01), std::max(_02, _10)), std::max(std::max(_12, _20), std::max(_21, _22)));
if (val >= vmax)//當前點和8鄰域相比是否為極大值
{
_00 = prevptr[c - step - 1]; _01 = prevptr[c - step]; _02 = prevptr[c - step + 1];
_10 = prevptr[c - 1]; _12 = prevptr[c + 1];
_20 = prevptr[c + step - 1]; _21 = prevptr[c + step]; _22 = prevptr[c + step + 1];
vmax = std::max(std::max(std::max(_00, _01), std::max(_02, _10)), std::max(std::max(_12, _20), std::max(_21, _22)));
if (val >= vmax)//當前點和上一層的8鄰域相比是否為極大值
{
_00 = nextptr[c - step - 1]; _01 = nextptr[c - step]; _02 = nextptr[c - step + 1];
_10 = nextptr[c - 1]; _12 = nextptr[c + 1];
_20 = nextptr[c + step - 1]; _21 = nextptr[c + step]; _22 = nextptr[c + step + 1];
vmax = std::max(std::max(std::max(_00, _01), std::max(_02, _10)), std::max(std::max(_12, _20), std::max(_21, _22)));
if (val >= vmax)//當前點和下一層的8鄰域相比是否為極大值
{
sift_wt _11p = prevptr[c], _11n = nextptr[c];
calculate = (val >= std::max(_11p, _11n));//當前點和上下兩層對應點相比是否為極大值
}
}
}
}
else { // val cant be zero here (first abs took care of zero), must be negative// val要不是正的，要不是負的，不可能是0
sift_wt vmin = std::min(std::min(std::min(_00, _01), std::min(_02, _10)), std::min(std::min(_12, _20), std::min(_21, _22)));
if (val <= vmin)
{
_00 = prevptr[c - step - 1]; _01 = prevptr[c - step]; _02 = prevptr[c - step + 1];
_10 = prevptr[c - 1]; _12 = prevptr[c + 1];
_20 = prevptr[c + step - 1]; _21 = prevptr[c + step]; _22 = prevptr[c + step + 1];
vmin = std::min(std::min(std::min(_00, _01), std::min(_02, _10)), std::min(std::min(_12, _20), std::min(_21, _22)));
if (val <= vmin)
{
_00 = nextptr[c - step - 1]; _01 = nextptr[c - step]; _02 = nextptr[c - step + 1];
_10 = nextptr[c - 1]; _12 = nextptr[c + 1];
_20 = nextptr[c + step - 1]; _21 = nextptr[c + step]; _22 = nextptr[c + step + 1];
vmin = std::min(std::min(std::min(_00, _01), std::min(_02, _10)), std::min(std::min(_12, _20), std::min(_21, _22)));
if (val <= vmin)
{
sift_wt _11p = prevptr[c], _11n = nextptr[c];
calculate = (val <= std::min(_11p, _11n));
}
}
}
}
...//後面還有求亞畫素精度極值，梯度方向直方圖
}
1.2.2 改進關鍵點位置的精度
一個連續函式被取樣時，它真正的最大值或最小值實際上可能位於樣本點之間。得到接近真實極值（亞畫素精度）的一種方法時，首先再數字函式中的每個極值點處擬合一個內插函式，然後在內插後的函式中查詢改進精度後的極值位置。SIFT用D(x,y,σ)的泰勒級數展開的線性項和二次項，把原點移至被檢測的樣本點。這個公式的向量形式為
式中，D及其導數是在這個樣本點處計算的，x=(x,y,σ)T是這個樣本的偏移量，▽是我們熟悉的梯度運算元，
H是海森矩陣，
取D(x)關於x的導數並令其為零，可求得極值的位置x，即
就像在部落格1.2節介紹的那樣，海森矩陣和D的梯度是用相鄰畫素點的差來近似的，f'(x)=(f(x+1)-f(x-1))/2，二階導數f''(x)=f(x+1)+f(x-1)-2f(x)。得到的3×3線性方程組加u三上很容易求解。如果偏移量x在其任意維度上都大於0.5，那麼可以得出結論：極值靠近另一個樣本點，在這種情況下，改變樣本點，並對改變後的樣本點進行內插。最後的偏移量x被新增到其樣本點的位置，得到極值的內插後的估計位置。
SIFT使用極值位置函式值D(x)來剔除具有低對比度的不穩定極值，其中D(x)為
Lowe的實驗結果稱，若所有影像的值都在區間[0,1]內，則D(x)小於0.03的任何極值都會被剔除。這剔除了具有低對比度和/或區域性化較差的關鍵點。
1.2.3 消除邊緣響應
回顧邊緣檢測的部落格可知，使用高斯差會得到影像中的邊緣。但SIFT中我們感興趣的關鍵點是“角狀”特徵，這些特徵更加區域性化。因此，消除了邊緣導致的灰度過渡。為了量化邊和角之間的差，我們研究區域性曲折度。邊在一個方向上由高曲折度表徵，在正交方向上由低曲折度表徵。影像中某點的曲折度可由該點處的一個2×2海森矩陣算出。因此，要在標量空間中的任何一層計算DoG的區域性曲折度，可在該層中計算D的海森矩陣：
H的特徵值與D的曲折度成正比。如對哈里斯-斯蒂芬斯(HS)角檢測器的解釋那樣，我們可以避免直接計算特徵值，H的跡等於特徵值之和，H的行列式等於特徵值之積。設α和β分別是H的最大特徵值和最小特徵值，有
注意，H是對稱的，且大小為2×2。若行列式是負的，則不同曲折度具有不同的符號，且討論的關鍵點不可能是一個極值，因此要丟棄它。令r表示最大特徵值和最小特徵值之比，即α=rβ，則
當特徵值相等時，上式出現最小值，當r大於1時，隨著r的增大而增大。因此，要檢查低於某個閾值的r的主曲折度之比，只需檢查
這個計算很簡單。Lowe報告的實驗結果中使用了值r=10，這意味著消除了曲折度之比大於10的關鍵點。
OpenCV原始碼adjustLocalExtrema
包括計算亞畫素精度極值+剔除低對比度的不穩定極值+剔除曲折度較大的點，基於Lowe論文的Section 4
static bool adjustLocalExtrema(
const std::vector<Mat>& dog_pyr, KeyPoint& kpt, int octv,
int& layer, int& r, int& c, int nOctaveLayers,
float contrastThreshold, float edgeThreshold, float sigma
)
{
CV_TRACE_FUNCTION();
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 = 0;
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];
//D(x,y,σ)在樣本點點的導數，用中心差分計算
Vec3f dD((img.at<sift_wt>(r, c + 1) - img.at<sift_wt>(r, c - 1)) * deriv_scale,
(img.at<sift_wt>(r + 1, c) - img.at<sift_wt>(r - 1, c)) * deriv_scale,
(next.at<sift_wt>(r, c) - prev.at<sift_wt>(r, c)) * deriv_scale);
// 基於中心差分的二階導數
float v2 = (float)img.at<sift_wt>(r, c) * 2;
float dxx = (img.at<sift_wt>(r, c + 1) + img.at<sift_wt>(r, c - 1) - v2) * second_deriv_scale;
float dyy = (img.at<sift_wt>(r + 1, c) + img.at<sift_wt>(r - 1, c) - v2) * second_deriv_scale;
float dss = (next.at<sift_wt>(r, c) + prev.at<sift_wt>(r, c) - v2) * second_deriv_scale;
float dxy = (img.at<sift_wt>(r + 1, c + 1) - img.at<sift_wt>(r + 1, c - 1) -
img.at<sift_wt>(r - 1, c + 1) + img.at<sift_wt>(r - 1, c - 1)) * cross_deriv_scale;
float dxs = (next.at<sift_wt>(r, c + 1) - next.at<sift_wt>(r, c - 1) -
prev.at<sift_wt>(r, c + 1) + prev.at<sift_wt>(r, c - 1)) * cross_deriv_scale;
float dys = (next.at<sift_wt>(r + 1, c) - next.at<sift_wt>(r - 1, c) -
prev.at<sift_wt>(r + 1, c) + prev.at<sift_wt>(r - 1, c)) * cross_deriv_scale;
//海森矩陣
Matx33f H(dxx, dxy, dxs,
dxy, dyy, dys,
dxs, dys, dss);
//取D(X)關於X的導數並令其為零，可求得極值的位置X, X = (x, y, σ)T
//即，HX=-dD
Vec3f X = H.solve(dD, DECOMP_LU);
xi = -X[2];
xr = -X[1];
xc = -X[0];
//如果由泰勒級數插值得到的三個座標的偏移量都小於0.5，說明已經找到特徵點，則退出迭代
if (std::abs(xi) < 0.5f && std::abs(xr) < 0.5f && std::abs(xc) < 0.5f)
break;
//如果三個座標偏移量中任意一個大於一個很大的數，則說明該極值點不是特徵點，函式返回
if (std::abs(xi) > (float)(INT_MAX / 3) ||
std::abs(xr) > (float)(INT_MAX / 3) ||
std::abs(xc) > (float)(INT_MAX / 3))
return false;
//由上面得到的偏移量重新定義插值中心的座標位置,用於下次迭代，上面偏移量小於0.5的話，四捨五入還在layer層的(r，c)點，就不用下面的幾句
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 //如果迭代次數超過最大迭代次數，這個極值點不是特徵點
if (i >= SIFT_MAX_INTERP_STEPS)
return false;
// 剔除具有低對比度的不穩定極值 D(X)=D+0.5*dD*X
{
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<sift_wt>(r, c + 1) - img.at<sift_wt>(r, c - 1)) * deriv_scale,
(img.at<sift_wt>(r + 1, c) - img.at<sift_wt>(r - 1, c)) * deriv_scale,
(next.at<sift_wt>(r, c) - prev.at<sift_wt>(r, c)) * deriv_scale);
float t = dD.dot(Matx31f(xc, xr, xi));
contr = img.at<sift_wt>(r, c) * img_scale + t * 0.5f;
if (std::abs(contr) * nOctaveLayers < contrastThreshold)//當響應值D(X)小於contrastThreshold時剔除
return false;
// 剔除曲折度較大的點，利用Hessian的跡和行列式計算曲折度，和上面不同，下面的Hessian是2×2的
// principal curvatures are computed using the trace and det of Hessian
float v2 = img.at<sift_wt>(r, c) * 2.f;
float dxx = (img.at<sift_wt>(r, c + 1) + img.at<sift_wt>(r, c - 1) - v2) * second_deriv_scale;
float dyy = (img.at<sift_wt>(r + 1, c) + img.at<sift_wt>(r - 1, c) - v2) * second_deriv_scale;
float dxy = (img.at<sift_wt>(r + 1, c + 1) - img.at<sift_wt>(r + 1, c - 1) -
img.at<sift_wt>(r - 1, c + 1) + img.at<sift_wt>(r - 1, c - 1)) * cross_deriv_scale;
float tr = dxx + dyy;
float det = dxx * dyy - dxy * dxy;
//若行列式是負的，則不同曲折度具有不同的符號，且討論的關鍵點不可能是一個極值，因此要丟棄它。曲折度大就說明兩個特徵值相差大，就是邊緣，丟棄
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;
kpt.response = std::abs(contr);
return true;
}
1.3 計算關鍵點方向
在之前過程中，我們得到了SIFT認為穩定的關鍵點。因為我們知道每個關鍵點在尺度空間中的位置，因此實現了尺度獨立性。下一步是根據區域性性質為每個關鍵點分配方向，實現影像旋轉的不變性。對此，SIFT使用了一種簡單的方法：
- 使用關鍵點的尺度來選擇最接近該尺度的高斯平滑影像L，然後採集關鍵點所在高斯金字塔影像3σ鄰域視窗內畫素的梯度幅度M(x,y)和方向角θ(x,y):
- 方向直方圖由每個關鍵點鄰域的樣本點的梯度方向形成。直方圖有36個容器，覆蓋了影像平面上360°的方向範圍（或者每45度一個柱共8個柱)。新增到直方圖中的每個樣本按其梯度幅度和一個圓形高斯函式加權（也就是說灰度差大的方向對關鍵點方向的影響更大，也增強特徵點近的鄰域點對關鍵點方向的作用，並減少突變的影響），圓形高斯函式的標準差是關鍵點尺度的1.5倍。
- 為了得到更精確的方向，通常還可以對離散的梯度直方圖進行插值擬合。具體而言，關鍵點的方向可以由和主峰值最近的三個柱值透過拋物線插值得到。
- 在梯度直方圖中，當存在一個相當於主峰值80%能量的柱值時，則可以在相同的位置以相同的尺度建立另一個關鍵點，關鍵點的方向不同。SIFT僅為15%左右的多方向點分配多個方向，但它們對影像匹配的貢獻非常明顯。
OpenCV原始碼calcOrientationHist
// Computes a gradient orientation histogram at a specified pixel
static float calcOrientationHist(
const Mat& img, Point pt, int radius,
float sigma, float* hist, int n
)
{
CV_TRACE_FUNCTION();
int i, j, k, len = (radius*2+1)*(radius*2+1);
float expf_scale = -1.f/(2.f * sigma * sigma);
cv::utils::BufferArea area;
float *X = 0, *Y = 0, *Mag, *Ori = 0, *W = 0, *temphist = 0;
area.allocate(X, len, CV_SIMD_WIDTH);
area.allocate(Y, len, CV_SIMD_WIDTH);
area.allocate(Ori, len, CV_SIMD_WIDTH);
area.allocate(W, len, CV_SIMD_WIDTH);
area.allocate(temphist, n+4, CV_SIMD_WIDTH);
area.commit();
temphist += 2;
Mag = X;
for( i = 0; i < n; i++ ) //初始化直方圖
temphist[i] = 0.f;
for( i = -radius, k = 0; i <= radius; i++ )//兩個for迴圈，計算(radius*2+1)*(radius*2+1)個dx，dy和w
{
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<sift_wt>(y, x+1) - img.at<sift_wt>(y, x-1));//計算相對值，不用除以2了
float dy = (float)(img.at<sift_wt>(y-1, x) - img.at<sift_wt>(y+1, x));
X[k] = dx; Y[k] = dy; W[k] = (i*i + j*j)*expf_scale;//這裡的W為高斯函式的e的指數
k++;
}
}
len = k;
// 計算鄰域中所有元素的高斯加權值W，梯度幅角Ori 和梯度幅值Mag
// compute gradient values, orientations and the weights over the pixel neighborhood
cv::hal::exp32f(W, W, len);
cv::hal::fastAtan2(Y, X, Ori, len, true);
cv::hal::magnitude32f(X, Y, Mag, len);
k = 0;
for( ; k < len; k++ )//漲知識了，我每次都要寫老長的判斷，判斷是在第幾象限
{
int bin = cvRound((n/360.f)*Ori[k]);//判斷鄰域畫素的梯度幅角屬於n個柱體的哪一個，和Lowe推薦的方法不同，它選了最近的柱體，然後再平滑
//如果超出範圍，則利用圓周迴圈確定其真正屬於的那個柱體，比如dy=-1，dx=1，Ori=-45°，假設n=8，bin=-1，bin+=n，就是第7個柱體
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];
// 平滑直方圖
i = 0;
for( ; 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;
}
1.4 關鍵點描述
經過前面的步驟，已經為關鍵點分配了影像位置、尺度和方向，進而為這三個變數提供了不變性。下一步是為圍繞每個明顯不同的關鍵點的區域性區域計算一個描述子，同時對尺度、方向、光照和影像視點的變化是儘可能不變的。基本思想是能用這些描述子來識別兩幅或多幅影像中區域性區域的匹配（相似性）。特徵描述子的生成大致有三個步驟：
1. 校正旋轉主方向，確保旋轉不變性。為了保證特徵向量的旋轉不變性，要以特徵點為中心，在附近鄰域內將座標軸旋轉θθ（特徵點的主方向）角度，即將座標軸旋轉為特徵點的主方向。旋轉後鄰域內畫素的新座標為：
2. 生成描述子，最終形成一個128維的特徵向量。
以關鍵點為中心取16×16畫素的一個區域，並使用畫素差計算這個區域中每個點處的梯度幅度和方向。然後使用標準差等於這個區域大小一般的高斯加權函式來分配一個權重，將這個權重乘以每個點處的梯度幅度。高斯加權函式在圖中顯示為一個圓，權重隨著到中心的距離增大而減小。這個函式的目的是，減少位置變化很小時描述子的突然變化。
將16×16的區域分為16個4×4區域，每個4×4區域中有16個方向，將所有的梯度方向量化為8個相隔45°的方向。SIFT不是把一個方向值分配給其最接近的容器，而是進行內插運算，根據這個值到每個容器中心的距離，按比例地在所有容器中分配直方圖的輸入。例如，容器的中心為[22.5, 67.5, 112.5, ..., 337.5]，每個輸入需要乘以權重1-d，d是從這個值到容器中心的最短距離，度量單位是直方圖間隔，最大的可能距離是1。假設某個方向值是22.5°，到第一個容器中心的距離是0，因此將滿輸入分配給這個容器。到下一個容器的距離是1，1-d=0，所以不分配給下一個容器，所有容器都是如此。假設某個方向值是45°，分配給第一個容器1/2，分配給第二個容器1/2。採用這種方法，為每個容器得到一個技術的比例分數，從而避免“邊界”效應，即描述子方向的細微變化導致的為不同容器分配到同一個值的效應。
16個4×4區域得到16個直方圖，將直方圖的8個方向顯示為一個小向量簇，其中每個向量的長度等於其對應容器的值。所以，一個描述子由4×4陣列組成，每個陣列包含8個方向值。在SIFT中，這描述子資料被組織為一個128維的向量。
3. 歸一化處理，將特徵向量長度進行歸一化處理，進一步去除光照的影響。
為了降低光照的影響，對特徵向量進行了兩階段的歸一化處理。首先，透過將每個分量除以向量範數，把向量歸一化為單位長度。由每個畫素值乘以一個常數所引起的影像對比度的變化，將以相同的常數乘以梯度，因此對比度的變化將被第一次歸一化抵消。每個畫素加上一個常量導致的亮度變化不會影響梯度值，因為梯度值是根據畫素插值計算的。因此描述子對光照的仿射變換是不變的。然後，也可能出現攝像機飽和等導致的非線性光照變化。這類變化會導致某些梯度的相對幅值的較大變化，但它們幾乎不會影響梯度方向。SIFT透過對歸一化特徵向量進行閾值處理，降低了較大梯度值的影響，使所有分量都小於實驗確定的值0.2。閾值處理後，特徵向量被重新歸一化為單位向量。
2 opencv中的SIFT
2.1 SIFT使用示例
#include <opencv2/opencv.hpp>
using namespace cv;
using namespace std;
int main() {
Mat src1 = imread("./1.jpg", 0);
Mat src2 = imread("./2.jpg", 0);
Mat img1 = imread("./2.jpg");
Mat img2 = imread("./2.jpg");
SIFT S;
Ptr<SIFT> pSIFT = SIFT::create(55);
vector<KeyPoint> points1, points2;
pSIFT->detect(src1, points1);
pSIFT->detect(src2, points2);
drawKeypoints(src1, points1, img1, Scalar(0,255,0));
drawKeypoints(src2, points2, img2, Scalar(0, 255, 0));
Mat discriptions1, discriptions2;
pSIFT->compute(src1, points1, discriptions1);
pSIFT->compute(src2, points2, discriptions2);
Ptr<DescriptorMatcher> matcher = DescriptorMatcher::create("BruteForce");
vector<DMatch> matches;
matcher->match(discriptions1, discriptions2, matches);
Mat match_img;
drawMatches(src1, points1, src2, points2, matches, match_img);
imshow("markImg", match_img);
waitKey(0);
return 0;
}
2.2 原始碼分析
SIFT建構函式，查詢使用者是否提供關鍵點，然後呼叫以下函式
- 建立初始影像SIFT::createInitialImage()
- 建立高斯金字塔SIFT::buildGaussianPyramid()√
- 建立DoG金字塔 SIFT::buildDoGPyramid()√
- 在DoG中查詢尺度空間極值點SIFT::findScaleSpaceExtrema()√，會呼叫下面兩個函式
（1）對找到的極值點進行曲線插值擬合，並過濾 adjustLocalExtrema()√
（2）計算方向直方圖 calcOrientationHist()√
- 計算描述子 calcDescriptors()
（1）計算sift描述子calcSIFTDescriptor()
上面在講原理的時候已經寫了一些原始碼註釋（打√），下面是剩下的註釋。
2.2.1 SIFT::createInitialImage()
其中σinit是第0層的尺度，σpre是被相機模糊後的尺度。在Lowe的論文使用了σ0=1.6，S=3。
為了儘可能多的保留原始影像資訊，一般需要對原始影像進行擴大兩倍取樣，即升取樣，從而生成一組取樣圖octave_1,此組取樣圖的第一層的模糊引數為：
也就是說，輸入的影像已經是被模糊了（0.5），但你想要更模糊一點（1.6）的作為base。如果要先擴大影像尺寸的話，就先升取樣後再用1.25的σ進行卷積。如果不升取樣的話，就用sqrt(1.6*1.6-0.5*0.5)對原圖做卷積，然後將結果作為base。
static Mat createInitialImage( const Mat& img, bool doubleImageSize, float sigma )
{
Mat gray, gray_fpt;
if( img.channels() == 3 || img.channels() == 4 )
{ //如果輸入影像是彩色影像，則需要轉換成灰度影像
cvtColor(img, gray, COLOR_BGR2GRAY);
gray.convertTo(gray_fpt, DataType<sift_wt>::type, SIFT_FIXPT_SCALE, 0);//調整影像的畫素資料型別
}
else
img.convertTo(gray_fpt, DataType<sift_wt>::type, SIFT_FIXPT_SCALE, 0);
float sig_diff;
if( doubleImageSize )//如果需要擴大影像的長寬尺寸
{
//SIFT_INIT_SIGMA為0.5，即輸入影像的尺度，SIFT_INIT_SIGMA×2 = 1.0，即影像擴大2倍以後的尺度
sig_diff = sqrtf( std::max(sigma * sigma - SIFT_INIT_SIGMA * SIFT_INIT_SIGMA * 4, 0.01f) );
Mat dbl;
#if DoG_TYPE_SHORT
resize(gray_fpt, dbl, Size(gray_fpt.cols*2, gray_fpt.rows*2), 0, 0, INTER_LINEAR_EXACT);
#else
resize(gray_fpt, dbl, Size(gray_fpt.cols*2, gray_fpt.rows*2), 0, 0, INTER_LINEAR);//利用雙線性插值法把影像的長寬都擴大2倍
#endif
Mat result;
GaussianBlur(dbl, result, Size(), sig_diff, sig_diff);//對影像進行高斯平滑處理
return result;
}
else//如果不需要擴大影像的尺寸，sig_diff計算不同
{
sig_diff = sqrtf( std::max(sigma * sigma - SIFT_INIT_SIGMA * SIFT_INIT_SIGMA, 0.01f) );
Mat result;
GaussianBlur(gray_fpt, result, Size(), sig_diff, sig_diff);
return result;
}
}
1. 計算鄰域內每個畫素點的幅值和幅角
在前面的步驟中我們分別得到了關鍵點的位置（層，組，所在影像的橫縱座標），尺度，主方向資訊。所以我們定位關鍵點到其所屬的高斯金字塔中的那幅圖中。用關鍵點周圍(radius*2+1)*(radius*2+1)的鄰域（橘色矩形框區域）來描述這個關鍵點，即計算鄰域的幅值和幅角，其中
將這塊區域分成d×d個小區域，SIFT中d取4，每個小區域的寬度是3σ，這裡的σ是指相對於當前組第一層圖片來說的。
2. 校正旋轉主方向，確保旋轉不變性。
為了保證特徵向量的旋轉不變性，要以特徵點為中心，在附近鄰域內將座標軸旋轉θ（特徵點的主方向）角度，即將座標軸旋轉為特徵點的主方向。如下圖所示
旋轉後鄰域內畫素的新座標為：
3.建立三維直方圖
在計算特徵點描述符的時候，我們不需要精確知道鄰域內所有畫素的梯度幅值和幅角，我們只需要根據直方圖知道其統計值即可。如圖所示，三維直方圖是由d*d*n個長寬為1的單位立方體組成的，高對應鄰域畫素幅角的大小，把360度分成8等分。立方體的底就是特徵點的鄰域區域，該區域被劃分為4×4個子區域，鄰域中的畫素根據座標位置，把它們歸屬到這16個子區域中的一個，再根據該畫素幅角的大小，把它分到這8等分中的一份，這樣每個畫素點都能對應到其中的一個立方體內，三維直方圖建立起來了。
正方體的中心代表著該正方體，但是落入正方體內的鄰域畫素不可能都在中心，因此我們需要對上面的梯度幅值做進一步的處理，根據它對中心點位置的貢獻大小進行加權處理，即在正方體內，根據畫素點相對於正方體中心的距離，對梯度幅值做加權處理。所以，三維直方圖的值，也就是正方體的值需要下面四個步驟完成：
1）計算落入該正方體內的鄰域畫素的梯度幅值A
2）根據該畫素相對於特徵點的距離，對A進行高斯加權處理，得到B
3）根據該畫素相對於它所在的正方體的中心的貢獻大小，再對B進行加權處理，得到C
由於計算相對於正方體中心點的貢獻大小略顯繁瑣,例如畫素（0.25,0.25,0.25），中心點為（0.5,0.5,0.5），相對距離是（0.25,0.25,0.25），貢獻權重為(1-0.25)*(1-0.25)*(1-0.25)。因此在實際應用中，我們需要經過座標平移，把中心點平移到正方體的頂點上,這樣只要計算正方體內的點對正方體的8個頂點的貢獻即可。根據三線性插值法，對某個頂點的貢獻值是以該頂點和正方體內的點為對角線的兩個頂點，所構成的立方體的體積。也就是對8個頂點的貢獻分別為：
經過上述處理後我們得到128維描述子，在OpenCV中，一共對描述子進行了兩次歸一化處理。第一次使用以下公式對描述子進行歸一化處理，並對大於0.2的描述子進行截斷處理：
假設經過第一次歸一化後描述子用qi表示，第二次歸一化操作是
其中分母要大於FLT_EPSILON。
void calcSIFTDescriptor(
const Mat& img, Point2f ptf, float ori, float scl,
int d, int n, Mat& dstMat, int row
)
{
Point pt(cvRound(ptf.x), cvRound(ptf.y));//特徵點座標位置
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;// 1/45
float exp_scale = -1.f/(d * d * 0.5f); //高斯加權函式中的e指數的常數部分
float hist_width = SIFT_DESCR_SCL_FCTR * scl;//SIFT_DESCR_SCL_FCTR = 3.f，即3σ
int radius = cvRound(hist_width * 1.4142135623730951f * (d + 1) * 0.5f);//特徵點鄰域區域的半徑
// Clip the radius to the diagonal of the image to avoid autobuffer too large exception
radius = std::min(radius, (int)std::sqrt(((double) img.cols)*img.cols + ((double) img.rows)*img.rows));//避免鄰域過大
cos_t /= hist_width;//歸一化處理
sin_t /= hist_width;
//len為特徵點鄰域區域內畫素的數量，histlen為直方圖的數量，即特徵向量的長度，實際應為d×d×n，+2是為圓周迴圈留出一定的記憶體空間
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;
//開闢一段記憶體空間
cv::utils::BufferArea area;
float *X = 0, *Y = 0, *Mag, *Ori = 0, *W = 0, *RBin = 0, *CBin = 0, *hist = 0, *rawDst = 0;
area.allocate(X, len, CV_SIMD_WIDTH);
area.allocate(Y, len, CV_SIMD_WIDTH);
area.allocate(Ori, len, CV_SIMD_WIDTH);
area.allocate(W, len, CV_SIMD_WIDTH);
area.allocate(RBin, len, CV_SIMD_WIDTH);
area.allocate(CBin, len, CV_SIMD_WIDTH);
area.allocate(hist, histlen, CV_SIMD_WIDTH);
area.allocate(rawDst, len, CV_SIMD_WIDTH);
area.commit();
Mag = Y;
//直方圖陣列hist清零
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.;
}
//遍歷當前特徵點的鄰域範圍
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;
//把鄰域區域的原點從中心位置移到該區域的左下角，以便後面的使用。因為變數cos_t 和sin_t 都已
//進行了歸一化處理，所以原點位移時只需要加d/2即可。而再減0.5f的目的是進行座標平移，從而在
//三線性插值計算中，計算的是正方體內的點對正方體8個頂點的貢獻大小，而不是對正方體的中心點的
//貢獻大小。之所以沒有對角度obin 進行座標平移，是因為角度是連續的量，無需平移
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;//得到鄰域畫素點的位置座標
//確定鄰域畫素是否在d×d 的正方形內，以及是否超過了影像邊界
if( rbin > -1 && rbin < d && cbin > -1 && cbin < d &&
r > 0 && r < rows - 1 && c > 0 && c < cols - 1 )
{
//計算x和y方向的一階導數，這裡沒有除以2，因為沒有分母部分不影響後面所進行的歸一化處理
float dx = (float)(img.at<sift_wt>(r, c+1) - img.at<sift_wt>(r, c-1));
float dy = (float)(img.at<sift_wt>(r-1, c) - img.at<sift_wt>(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;//高斯加權函式中的e 指數部分
k++;//統計實際的鄰域畫素的數量
}
}
len = k;
cv::hal::fastAtan2(Y, X, Ori, len, true);//計算梯度幅角
cv::hal::magnitude32f(X, Y, Mag, len);//計算梯度幅值
cv::hal::exp32f(W, W, len);//計算高斯加權函式
k = 0;
for( ; k < len; k++ )
{
float rbin = RBin[k], cbin = CBin[k];//得到d×d 鄰域區域的座標，即三維直方圖的底內的位置
float obin = (Ori[k] - ori)*bins_per_rad;//得到幅角所屬的8等份中的某一個等份，即三維直方圖的高的位置
float mag = Mag[k]*W[k];//得到高斯加權以後的梯度幅值
//r0，c0和o0 為三維座標的整數部分，它表示屬於的哪個正方體
int r0 = cvFloor( rbin );//向下取整
int c0 = cvFloor( cbin );
int o0 = cvFloor( obin );
//rbin，cbin和obin 為三維座標的小數部分，即將中心點移動到正方體端點時C點在正方體內的座標
rbin -= r0;
cbin -= c0;
obin -= o0;
//如果角度o0 小於0 度或大於360 度，則根據圓周迴圈，把該角度調整到0～360度之間
if( o0 < 0 )
o0 += n;
if( o0 >= n )
o0 -= n;
// histogram update using tri-linear interpolation
//根據三線性插值法，計算該畫素對正方體的8個頂點的貢獻大小，即8個立方體的體積，這裡還需要乘以高斯加權後的梯度值mag
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;
//8個頂點對應於座標平移前的8個直方圖的正方體，對其進行累加求和
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
//由於圓周迴圈的特性，縱向的n+2個正方體，使得第n個和第n+1個的值加到第0個和第1個上，而底部的（d+2）*（d+2）則是捨棄多出來的2
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++ )
rawDst[(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;
k = 0;
for( ; k < len; k++ )
nrm2 += rawDst[k]*rawDst[k];//將所有特徵向量的平方和累加起來
float thr = std::sqrt(nrm2)*SIFT_DESCR_MAG_THR;//累加平方和*0.2為閾值
i = 0, nrm2 = 0;
for( ; i < len; i++ )
{
float val = std::min(rawDst[i], thr);//大於閾值的值用thr代替
rawDst[i] = val;
nrm2 += val*val;
}
nrm2 = SIFT_INT_DESCR_FCTR/std::max(std::sqrt(nrm2), FLT_EPSILON);//截斷後的用於歸一化的值
k = 0;
if( dstMat.type() == CV_32F )
{
float* dst = dstMat.ptr<float>(row);
for( ; k < len; k++ )
{
dst[k] = saturate_cast<uchar>(rawDst[k]*nrm2);//歸一化
}
}
else // CV_8U
{
uint8_t* dst = dstMat.ptr<uint8_t>(row);
for( ; k < len; k++ )
{
dst[k] = saturate_cast<uchar>(rawDst[k]*nrm2);
}
}
}
參考
1. 岡薩雷斯《數字影像處理（第四版）》Chapter 11（所有圖片可在連結中下載）
2. 影像處理基礎（六）基於 LOG 的 blob 興趣點檢測
4. 影像處理基礎 （四）邊緣提取之 LOG 和 DOG 運算元
6. sift程式詳解